---
title: "ICC final version"
author: "Anne Christine Bischops"
output: html_document
---

# Data cleaning and merging
Load libraries.
```{r load packages}

library(tidyverse) 
library(dplyr)# already included in tidyverse, but can alternatively be installed on its own
library(ggrepel) # to avoid text labels in ggplot from overlapping
library(srvyr)  # survey package that also works with dplyr 
library(readstata13) # foreign only reads Stata 12 or lower
library(tableone) # Creates a table 1 (summary characteristics)
library(lme4)
library(spatstat)

memory.limit(size=50000) #increases the limit in force on the total allocation due to the large dataset

```


The code chunk below takes the DLHS and AHS datasets that were cleaned in Stata and contain sampling weights, cleans them further, drops irrelevant columns, and then appends DLHS and AHS to yield a dataset for analysis ("india.cvd"). The code also adds age weights from the WHO standard population, which I calculated separately for adults in an excel file. These weights are multiplied by the sampling weights (sweight) to yield a weight that can be used for age standardization (sworld_weight). 

```{r Cleaning and merging data, eval=FALSE}
# 1. Import DLHS data which was previously cleaned in Stata (detailed documentation in the repository) ########################################################################################
dlhs <- read.dta13("DLHS_w3_v2.dta")
dlhs <- as_tibble(dlhs)

# 2. Drop irrelevant columns in DLHS ##############################################################################################################
dlhs <- dplyr::select(dlhs, state, state_dist, female, educ, age1, rural, q_ai_nat_rurb, q_ai_nat,
                      caste, hindu, muslim, christian, sikh, buddhist, jain, noreligion, otherreligion,
                      dropage, pregnant, 
                      married,
                      
                      BMI,
                      bpsyst_avg, bpsyst_frst, bpsyst_sec, htn_treated, htn_broad_avg, bpdiast_frst, bpdiast_sec,
                      glucose, diab_broad, fasted, 
                      hv98,  # this is the smoking var
                      
                      hhid, DLHS, sweight, psu, # vars needed for weighting
                      ai_NAT, ai_NAT_rurb, htn_narrow_avg,diab_narrow,diab_narrow_adj, diab_narrow_adj_nf, diab_narrow, diab_narrow_nf,adulthh) 

# 3. Filter out those <18 or pregnant  ##########################################################################################################
dlhs <- dplyr::filter(dlhs, dropage == 0) # only those >18 and with non-missing age
dlhs <- dplyr::filter(dlhs, pregnant == 0)  # only keep those who aren't pregnant (they didn't measure glucose in pregnant women anyway); DLHS has no missings in the pregnant variable
dlhs <- dplyr::select(dlhs, -dropage, -pregnant)
dlhs.cvd <- dlhs
dlhs.cvd <- dlhs.cvd %>% 
  dplyr::rename(age = age1, 
                smoke = hv98) %>% 
  mutate(hhid=as.character(hhid))


# 4. Import AHS data which was previously cleaned in Stata (detailed documentation in the repository)##################################################################################
ahs <- read.dta13("AHS_w3_v2.dta")
ahs <- as_tibble(ahs)

# 5. Drop irrelevant columns in AHS ######################################################################################################
ahs <- dplyr::select(ahs, state, state_dist, female, educ, age1, rural, q_ai_nat_rurb, q_ai_nat,
                      caste, hindu, muslim, christian, sikh, buddhist, jain, noreligion, otherreligion,
                      dropage, pregnant,
                      married,adulthh,
                      
                      BMI,
                      bpsyst_avg, bpsyst_frst, bpsyst_sec, htn_treated, htn_broad_avg, bpdiast_frst, bpdiast_sec,
                      glucose, diab_broad, 
                      smoke,  # this is the equivalent to hv98 in dlhs
                      
                      fid, DLHS, sweight, psu,  # vars needed for weighting; fid is the hhid equivalent
                      ai_NAT, ai_NAT_rurb, htn_narrow_avg, diab_narrow, diab_narrow_adj, diab_narrow_adj_nf,diab_narrow_nf) 


ahs.cvd <- ahs %>% 
  dplyr::rename(age = age1) %>% 
  mutate(hhid=as.character(fid))  # fid is the household ID in the AHS. 


# 6. Filter out those pregnant and <18 ###################################################################################################

ahs.cvd <- dplyr::filter(ahs.cvd, dropage == 0)
ahs.cvd <- dplyr::filter(ahs.cvd, pregnant == 0 | is.na(pregnant)==TRUE)  
ahs.cvd <- dplyr::select(ahs.cvd, -dropage, -pregnant)

# 7. Merge AHS AND DLHS  ##############################################################################
dlhs.cvd$state <- as.factor(dlhs.cvd$state)
ahs.cvd$state <- as.factor(ahs.cvd$state)
#smoke variable needs to be converted to numeric in ahs
ahs.cvd$smoke <-as.numeric(ahs.cvd$smoke)
india.cvd <- bind_rows(dlhs.cvd, ahs.cvd) # ahs.cvd has two less columns than dlhs.cvd because the fasted variable and psu is only in DLHS. R sets this to missing for AHS which is fine. 


# 8. Convert dbl to factors & clean  #################################################################
india.cvd <- dplyr::select(india.cvd, glucose, BMI, age, bpsyst_avg, bpsyst_frst, bpsyst_sec, bpdiast_frst, bpdiast_sec, sweight, ai_NAT, ai_NAT_rurb, DLHS, everything())

india.cvd[,12:36] <- lapply(india.cvd[,12:36], factor)  # factors all vars but the first 10 

india.cvd <- dplyr::mutate(india.cvd, educ=fct_recode(educ,
  "<Primary" = "1", 
  "Primary" = "2", 
  "Middle" = "3",
  "Secondary" = "4",
  "High School" = "5",
  ">High School" = "6"))

# This step is necessary because the DLHS and AHS states came without labels (just IDs)
india.cvd <- mutate(india.cvd,
  state = ifelse(state==2, "Himachal Pradesh", #ex_state_ind previously, error, only MP
  ifelse(state==3, "Punjab",
  ifelse(state==4, "Chandigarh",
  ifelse(state==6, "Haryana",
  ifelse(state==7, "Delhi",
  ifelse(state==11, "Sikkim",
  ifelse(state==12, "Arunachal Pradesh",
  ifelse(state==13, "Nagaland",
  ifelse(state==14, "Manipur",
  ifelse(state==15, "Mizoram",
  ifelse(state==16, "Tripura",
  ifelse(state==17, "Meghalaya",
  ifelse(state==19, "West Bengal",
  ifelse(state==25, "Daman and Diu",
  ifelse(state==27, "Maharashtra",
  ifelse(state==28, "Andhra Pradesh",
  ifelse(state==29, "Karnataka",
  ifelse(state==30, "Goa",
  ifelse(state==32, "Kerala",
	ifelse(state==33, "Tamil Nadu",
	ifelse(state==34, "Puducherry",
	ifelse(state==35, "Andaman and Nicobar",
	ifelse(state==36, "Telangana",
	ifelse(state==5, "Uttarakhand",
	ifelse(state==8, "Rajasthan",
	ifelse(state==9, "Uttar Pradesh",
	ifelse(state==10, "Bihar",
	ifelse(state==18, "Assam",
	ifelse(state==20, "Jharkhand",
	ifelse(state==21, "Odisha",
	ifelse(state==22, "Chhattisgarh", "Madhya Pradesh"))))))))))))))))))))))))))))))))
india.cvd$state <- as.factor(india.cvd$state)

# 9. Create age and bmi group ##############################################################################
india.cvd <- dplyr::mutate(india.cvd, age_grp = ifelse(age<=25, "18-25", 
                                                         ifelse(age>25 & age<=35, "26-35",
                                                                ifelse(age>35 & age<=45, "36-45",
                                                                       ifelse(age>45 & age<=55, "46-55",
                                                                              ifelse(age>55 & age<=65, "56-65", ">65")))))) 

india.cvd$age_grp <- factor(india.cvd$age_grp, levels = c("18-25", "26-35", "36-45", "46-55", "56-65", ">65"))
india.cvd <- within(india.cvd, age_grp <- relevel(age_grp, ref = "18-25"))

#create 5 year age groups
india.cvd <- dplyr::mutate(india.cvd, age_grp2 = ifelse(age>=15 &age<=19 , "15-19", 
                                                         ifelse(age>19 & age<=24, "20-24",
                                                                ifelse(age>24 & age<=29, "25-29",
                                                                       ifelse(age>29 & age<=34, "30-34",
                                                                              ifelse(age>34 & age<=39, "35-39", 
                                                                                     ifelse(age>39 & age<=44, "40-44",
                                                                                            ifelse(age>44 & age<=49, "45-49",
                                                                                                   ifelse(age>49 & age<=54, "50-54",
                                                                                                          ifelse(age>54 & age<=59, "55-59",
                                                                                                                 ifelse(age>59 & age<=64, "60-64",
                                                                                                                        ifelse(age>64 & age<=69, "65-69",
                                                                                                                               ifelse(age>69 & age<=74, "70-74",
                                                                                                                                      ifelse(age>74 & age<=79, "75-79","80+")))))))))))))) 

#create numerized age_grp5
india.cvd <- dplyr::mutate(india.cvd, age_5yr_2 = ifelse(age>=15 &age<=19 , 1, 
                                                         ifelse(age>19 & age<=24, 2,
                                                                ifelse(age>24 & age<=29, 3,
                                                                       ifelse(age>29 & age<=34, 4,
                                                                              ifelse(age>34 & age<=39, 5, 
                                                                                     ifelse(age>39 & age<=44, 6,
                                                                                            ifelse(age>44 & age<=49, 7,
                                                                                                   ifelse(age>49 & age<=54, 8,
                                                                                                          ifelse(age>54 & age<=59, 9,
                                                                                                                 ifelse(age>59 & age<=64, 10,
                                                                                                                        ifelse(age>64 & age<=69, 11,
                                                                                                                               ifelse(age>69 & age<=74, 12,
                                                                                                                                      ifelse(age>74 & age<=79, 13,14)))))))))))))) 

#bmi variables ICC
india.cvd$bmi<-india.cvd$BMI
india.cvd <- india.cvd %>% 
  mutate(bmicat2=ifelse(bmi<16,1,
                          ifelse(bmi>=16& bmi<18.5,2,
                          ifelse(bmi>=18.5& bmi<25,3,
                                              ifelse(bmi>=25 &bmi<30,4,
                                                     ifelse(bmi>=30,5,NA)))))) 

#bmi categories according to "Appropriate body-mass index for Asian populations and its implications for policy and intervention strategies. Lancet 2004 "                                              
india.cvd <- india.cvd %>% 
  mutate(bmicat3=ifelse(bmi<16,0,
                        ifelse(bmi>=16 &bmi<18.5,1,
                          ifelse(bmi>=18.5& bmi<23,2,
                          ifelse(bmi>=23& bmi<25,3,
                                              ifelse(bmi>=25 &bmi<27.5,4,
                                                     ifelse(bmi>=27.5 & bmi<30,5,
                                                            ifelse(bmi>=30,6,NA))))))))

india.cvd$bmicat3 <-as.factor(india.cvd$bmicat3)

# Calculate obese
india.cvd <- dplyr::mutate(india.cvd, 
                            obese = ifelse(BMI>=30, "1", "0"))
india.cvd$obese <- as.factor(india.cvd$obese)

#create new overweight and obesity variables
india.cvd <- india.cvd %>%
  mutate(overweight23=ifelse(is.na(bmi)==T,NA,
                             ifelse(bmi>=23,1,0)),
        obese27=ifelse(is.na(bmi)==T,NA,
                     ifelse(bmi>=27.5,1,0)))

# 10. Create smoking variables
##########################################################################################################
# Current smoker:
india.cvd <- mutate(india.cvd, 
                   currsmoke = ifelse(smoke==1 | smoke==2 | smoke == "occational-smoker" | smoke == "usual-smoker", 1,
                                      ifelse(smoke==3 | smoke==4 | smoke == "ex-smoker" | smoke == "never-smoked", 0, NA)))
india.cvd$currsmoke <- as.factor(india.cvd$currsmoke)

# ever smoker
india.cvd <- mutate(india.cvd, 
                   eversmoke = ifelse(smoke==1 | smoke==2 | smoke==3 | smoke == "occational-smoker" | smoke == "usual-smoker" | smoke == "ex-smoker", 1,
                                      ifelse(smoke==4 | smoke == "never-smoked", 0, NA)))
india.cvd$eversmoke <- as.factor(india.cvd$eversmoke)


# 11. Clean up caste, relig, create male, urban, and labelled vars ######################################
# AHS contains only ST, SC, and other, while DLHS has ST, SC, OBC, and other. I'm merging OBC and other in DLHS so that it's comparable. 

india.cvd <- mutate(india.cvd, 
                     religion = ifelse(hindu=="yes", "Hindu", 
                                       ifelse(christian=="yes", "Christian",
                                              ifelse(muslim=="yes", "Muslim",
                                                     ifelse(sikh=="yes", "Sikh", 
                                                            ifelse(buddhist=="yes", "Buddhist",
                                                                    ifelse(noreligion=="yes" | otherreligion=="yes" | jain=="yes", "Other", NA)))))))  
india.cvd$religion <- as.factor(india.cvd$religion)   

india.cvd <- mutate(india.cvd, 
                     caste = fct_recode(caste, 
                          "Scheduled caste" = "1",
                          "Scheduled tribe" = "2",
                          "Other"="3",
                          "Other"="4"))

india.cvd <- mutate(india.cvd, 
                     urban = ifelse(rural=="1", 0, 1),
                     male = ifelse(female=="1", 0, 1))
india.cvd$urban <- as.factor(india.cvd$urban)
india.cvd$male <- as.factor(india.cvd$male)

india.cvd <- india.cvd %>% 
  mutate(female_lab = as.factor(ifelse(is.na(female) == TRUE, NA, 
                                    ifelse(female == 1, "Female", "Male"))),
         urban_lab = as.factor(ifelse(is.na(rural) == TRUE, NA, 
                                    ifelse(rural == 1, "Rural", "Urban"))),
         q_ai_nat_rurb_lab = as.factor(ifelse(is.na(q_ai_nat_rurb) == TRUE, NA, 
                                    ifelse(q_ai_nat_rurb == 1, "Poorest", 
                                           ifelse(q_ai_nat_rurb == 5, "Richest", q_ai_nat_rurb)))))

india.cvd$q_ai_nat_rurb_lab <- as.factor(india.cvd$q_ai_nat_rurb_lab)
india.cvd <- within(india.cvd, q_ai_nat_rurb_lab <- relevel(q_ai_nat_rurb_lab, ref = "Poorest"))


# 12. Correct missings in diabetes and htn  ######################################################################################
india.cvd <- india.cvd %>% 
  mutate(diab_narrow = ifelse(is.na(glucose)==TRUE, NA, diab_narrow)) 
india.cvd$diab_narrow <- as.factor(india.cvd$diab_narrow)
india.cvd <- india.cvd %>% 
  mutate(htn_narrow_avg = ifelse(is.na(bpsyst_frst)==TRUE | is.na(bpsyst_sec)==TRUE | is.na(bpdiast_frst)==TRUE | is.na(bpdiast_sec)==TRUE, NA, htn_narrow_avg))
india.cvd$htn_narrow_avg <- as.factor(india.cvd$htn_narrow_avg)



# 13. Sampling weights  ######################################################################################
# Give missing sweight obs the average sweight (very few obs have a missing sweight)

india.cvd <- mutate(india.cvd, 
                     sweight = ifelse(is.na(sweight)==TRUE, mean(sweight, na.rm=TRUE), sweight))

#merge age standardization weight from GBD India pop into the dataset

agest_weight <- read.csv("GBDpopweights_2018-04-24-age_grp15-19.csv")
agest_weight$female <-as.factor(agest_weight$sex) #sex is coded as 1=female
agest_weight<- dplyr::select(agest_weight, -1)

india.cvd<- left_join(india.cvd, agest_weight) 
india.cvd<- mutate(india.cvd, 
                    sworld_weight_india = sweight*gbd_weight)

india.cvd <- mutate(india.cvd, 
                    sworld_weight_india = ifelse(is.na(gbd_weight)==TRUE, mean(sworld_weight_india, na.rm=TRUE), sworld_weight_india))


```


```{r create zones and dbl variables}

# Add dbl vars
india.cvd <- mutate(india.cvd, 
                    diab_narrow_dbl = as.numeric(diab_narrow)-1,
                    htn_treated_dbl = as.numeric(htn_treated)-1,
                    female_dbl = as.numeric(female)-1)



################## Zones as per: https://en.wikipedia.org/wiki/Administrative_divisions_of_India
india.cvd <- mutate(india.cvd, 
                    # Nothern
    zone = as.factor(ifelse(state=="Chandigarh" | state=="NCT of Delhi" | state=="Haryana" | state=="Himachal Pradesh" | state=="Punjab" | state=="Rajasthan", "North",
                    # Northeastern
           ifelse(state=="Assam" | state=="Arunachal Pradesh" | state=="Manipur" | state=="Meghalaya" | state=="Mizoram" | state=="Nagaland" | state=="Sikkim" | state=="Tripura", "Northeast",
                    # Central
           ifelse(state=="Chhattisgarh" | state=="Madhya Pradesh" | state=="Uttarakhand" | state== "Uttar Pradesh", "Central",
                    # Eastern
           ifelse(state=="Bihar" | state=="Jharkhand" | state=="Odisha" | state=="West Bengal", "East",
                    # Western
           ifelse(state=="Daman and Diu" | state=="Goa" | state=="Maharashtra", "West",
                    # Southern
           ifelse(state=="Andaman and Nicobar" | state=="Andhra Pradesh" | state=="Karnataka" | state=="Kerala" | state=="Puducherry" | state=="Tamil Nadu" | state=="Telangana", "South", NA))))))))


# create stratum ID
india.cvd$state_dist_str <- as.character(india.cvd$state_dist)
india.cvd$rural_str <- as.character(india.cvd$rural)
india.cvd <- mutate(india.cvd, 
                     stratumid = str_c(state_dist_str, rural_str, sep = "_")) 
india.cvd$stratumid <- as.factor(india.cvd$stratumid)

#create PSU ID:
india.cvd <- india.cvd %>% 
                   mutate(psu_str = as.character(psu), 
                          psuid = str_c(state_dist_str, psu_str, sep = "_"))
india.cvd$psuid <- as.factor(india.cvd$psuid)


# After trying all possible combinations, the below is the best way of specifying the survey design. # The point estimates are unaffected by any of the svy choices unless sweight is not included. The only impact is on the SEs and the SEs are very small in every possible scenario. I have decided not to include fpc as district is a stratum (rather than a sampling unit) and I'm fairly confident that usually <5% of PSUs in each districts were selected. The inclusion of stratumid does not affect the point estimates or CIs but I'm still included it in just in case. The fact that I can't include PSU-ID in the AHS states and the overall svydesign makes my CIs a bit too tight but I don't think it's worth worrying about (since CIs are very tight here anyway). 

```

ANALYSIS 

```{r Filter outcome variables}

india<-india.cvd #renamed (to have india.cvd as a backup version)

# restrict to non-missings for all outcomes (anemia,overweight,underweight,diab,htn)
#so no missings in BMI,blood glucose, BP, hemoglobine and sex
india<- filter(india, is.na(female)==F & is.na(bmi)==F & is.na(currsmoke)==F) #sex and BMI and csmoke
india <- filter(india, is.na(htn_narrow_avg)==F & is.na(diab_narrow)==F)

#1,103,476 obs

#change variable names as used in previous code versions
india$ex_htn_narrow_ind <-india$htn_narrow_avg
india$ex_diab_narrow_ind <- india$diab_narrow 
india$csmoke <-india$currsmoke
india$wealth_quintile_rurb<-india$q_ai_nat_rurb
india$marital <- india$married
india$age_10yr <-india$age_grp
india$d_id <-india$state_dist
india$asset_index <-india$ai_NAT_rurb
india$ex_state_ind <- india$state

```

```{r 1.To what degree do outcomes cluster in diffent social/geographic units}
#prep for ICC calculation

#regress outcome with random intercept for level and then extract variance components to calculate ICC

#write function for ICC calculation with anova method
calc.icc <- function(y) {
  vc <-VarCorr(y) 
  vc<-as.data.frame(vc) 
  (vc$vcov[1] / (vc$vcov[1] + vc$vcov[2]))} 

```

```{r Table 1, eval=FALSE}
#########Sample characteristics#######################################
#check whether all are factors (csmoke!)

table1names <- c("female", "ex_htn_narrow_ind","ex_diab_narrow_ind", "bmi","bmicat2","bmicat3","age","age_grp",
                "urban", "q_ai_nat_rurb","educ","married","csmoke")

totalmiss <- CreateTableOne(vars = table1names, data=india, includeNA=TRUE)
print(totalmiss)
totalnomiss <- CreateTableOne(vars = table1names, data=india, includeNA=F)
print(totalnomiss)

sexmiss <- CreateTableOne(vars = table1names, data=india, strata=c("female"), includeNA=T)
print(sexmiss)
sexnomiss <- CreateTableOne(vars = table1names, data=india, strata=c("female"), includeNA=F)
print(sexnomiss)

```

```{r 1.To what degree do outcomes cluster in diffent social/geographic units}
#I calculate an ICC for each outcome(diabetes,hypertension,smoking,overweight,obesity) for the HH/PSU/District/State level


#################Fit the random effects model and calculate the ICC###############
#response must be numeric
india$ex_diab_narrow_ind<-as.numeric(levels(india$ex_diab_narrow_ind))[india$ex_diab_narrow_ind]
india$ex_htn_narrow_ind<-as.numeric(levels(india$ex_htn_narrow_ind))[india$ex_htn_narrow_ind]
india$csmoke<-as.numeric(levels(india$csmoke))[india$csmoke]
#india$overweight23<-as.numeric(levels(india$overweight23))[india$overweight23]
#india$obese27<-as.numeric(levels(india$obese27))[india$obese27]


####table 2####################################################
#state level
diab_state <- lmer(ex_diab_narrow_ind ~ 1 + (1|state), data=india)
calc.icc(diab_state)

htn_state <- lmer(ex_htn_narrow_ind ~ 1 + (1|state), data=india)
calc.icc(htn_state)

csmoke_state <- lmer(csmoke ~ 1 + (1|state), data=india)
calc.icc(csmoke_state)

overweight_state <- lmer(overweight23 ~ 1 + (1|state), data=india)
calc.icc(overweight_state)

obese_state <- lmer(obese27 ~ 1 + (1|state), data=india)
calc.icc(obese_state)


#calculcate confidence intervals via bootstrapping

#Calculate the bootstrap distribution
boot_diab_state <- bootMer(diab_state, calc.icc, nsim=500)#minimum simulations 500
boot_htn_state <-bootMer(htn_state, calc.icc, nsim=500)
boot_csmoke_state <-bootMer(csmoke_state, calc.icc, nsim=500)
boot_overweight23_state <- bootMer(overweight_state, calc.icc, nsim=500)
boot_obese27_state <-bootMer(obese_state, calc.icc, nsim=500)

#Draw from the bootstrap distribution the usual 95% upper and lower confidence limits
quantile(boot_diab_state$t, c(0.025, 0.975))

quantile(boot_htn_state$t, c(0.025, 0.975))
quantile(boot_csmoke_state$t, c(0.025, 0.975))
quantile(boot_overweight23_state$t, c(0.025, 0.975))
quantile(boot_obese27_state$t, c(0.025, 0.975))

#district level###############################################################
diab_dist <- lmer(ex_diab_narrow_ind ~ 1 + (1|d_id), data=india)
calc.icc(diab_dist)

htn_dist <- lmer(ex_htn_narrow_ind ~ 1 + (1|d_id), data=india)
calc.icc(htn_dist)

csmoke_dist <- lmer(csmoke ~ 1 + (1|d_id), data=india)
calc.icc(csmoke_dist)

overweight_dist <- lmer(overweight23 ~ 1 + (1|d_id), data=india)
calc.icc(overweight_dist)

obese_dist<- lmer(obese27 ~ 1 + (1|d_id), data=india)
calc.icc(obese_dist)


#calculcate confidence intervals via bootstrapping
#Calculate the bootstrap distribution
memory.limit(size=50000)
boot_diab_dist <- bootMer(diab_dist, calc.icc, nsim=500) #minimum simulations 500
boot_htn_dist <-bootMer(htn_dist, calc.icc, nsim=500)
boot_csmoke_dist <-bootMer(csmoke_dist, calc.icc, nsim=500)
boot_overweight_dist <- bootMer(overweight_dist, calc.icc, nsim=500)
boot_obese_dist <-bootMer(obese_dist, calc.icc, nsim=500)

#Draw from the bootstrap distribution the usual 95% upper and lower confidence limits
quantile(boot_diab_dist$t, c(0.025, 0.975))

quantile(boot_htn_dist$t, c(0.025, 0.975))
quantile(boot_csmoke_dist$t, c(0.025, 0.975))
quantile(boot_overweight_dist$t, c(0.025, 0.975))
quantile(boot_obese_dist$t, c(0.025, 0.975))


##########################PSU level###############################################
diab_psu <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=india)
calc.icc(diab_psu)

htn_psu <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=india)
calc.icc(htn_psu)

csmoke_psu <- lmer(csmoke ~ 1 + (1|psuid), data=india)
calc.icc(csmoke_psu)

overweight_psu <- lmer(overweight23 ~ 1 + (1|psuid), data=india)
calc.icc(overweight_psu)


obese_psu<- lmer(obese27 ~ 1 + (1|psuid), data=india)
calc.icc(obese_psu)


#calculcate confidence intervals via bootstrapping
#Calculate the bootstrap distribution
boot_diab_psu <- bootMer(diab_psu, calc.icc, nsim=500) #minimum simulations 500
boot_htn_psu <-bootMer(htn_psu, calc.icc, nsim=500)
boot_csmoke_psu <-bootMer(csmoke_psu, calc.icc, nsim=500)
boot_overweight_psu <- bootMer(overweight_psu, calc.icc, nsim=500)
boot_obese_psu <-bootMer(obese_psu, calc.icc, nsim=500)

#Draw from the bootstrap distribution the usual 95% upper and lower confidence limits
quantile(boot_diab_psu$t, c(0.025, 0.975))

quantile(boot_htn_psu$t, c(0.025, 0.975))
quantile(boot_csmoke_psu$t, c(0.025, 0.975))
quantile(boot_overweight_psu$t, c(0.025, 0.975))
quantile(boot_obese_psu$t, c(0.025, 0.975))


####household level###################################################################
diab_hh <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=india)
calc.icc(diab_hh)

htn_hh <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=india) 
calc.icc(htn_hh)

csmoke_hh <- lmer(csmoke ~ 1 + (1|hh_id), data=india)
calc.icc(csmoke_hh)

overweight_hh <- lmer(overweight23 ~ 1 + (1|hh_id), data=india)
calc.icc(overweight_hh)

obese_hh<- lmer(obese27 ~ 1 + (1|hh_id), data=india)
calc.icc(obese_hh)


#calculcate confidence intervals via bootstrapping
#Calculate the bootstrap distribution
boot_diab_hh <- bootMer(diab_hh, calc.icc, nsim=500) #minimum simulations 500
boot_htn_hh <-bootMer(htn_hh, calc.icc, nsim=500)
boot_csmoke_hh <-bootMer(csmoke_hh, calc.icc, nsim=500)
boot_overweight_hh <- bootMer(overweight_hh, calc.icc, nsim=500)
boot_obese_hh <-bootMer(obese_hh, calc.icc, nsim=500)

#Draw from the bootstrap distribution the usual 95% upper and lower confidence limits
quantile(boot_diab_hh$t, c(0.025, 0.975))

quantile(boot_htn_hh$t, c(0.025, 0.975))
quantile(boot_csmoke_hh$t, c(0.025, 0.975))
quantile(boot_overweight_hh$t, c(0.025, 0.975))
quantile(boot_obese_hh$t, c(0.025, 0.975))


```

```{r figure 1: ICC calculation for each state and state map}

#Calculate household ICC for each state separately (for each outcome) and then create a map of India where you map state-level variation of household ICC

#####################################ICCs by state ######
#calculate ICCs for each state
#divide dataset into dataset by state

india$hh_id<-as.numeric(levels(india$hh_id))[india$hh_id]
#india$psuid<-as.numeric(levels(india$psuid))[india$psuid]

AndamanandNicobarIslands <- india %>%
  filter(state== "Andaman and Nicobar")
AndhraPradesh <- india %>%
  filter(state== "Andhra Pradesh")
ArunachalPradesh <- india %>%
  filter(state== "Arunachal Pradesh")
Assam <- india %>%
  filter(state== "Assam")
Bihar <- india %>%
  filter(state== "Bihar")
Chandigarh <- india %>%
  filter(state== "Chandigarh")
Chhattisgarh <- india %>%
  filter(state== "Chhattisgarh")
#DadraandNagarHaveli <- india %>%
  #filter(state== "Dadra and Nagar Haveli")
DamanandDiu <- india %>%
  filter(state== "Daman and Diu")
Delhi <- india %>%
  filter(state== "Delhi")
Goa <- india %>%
  filter(state== "Goa")
#Gujarat <- india %>%
  #filter(state== "Gujarat")
Haryana <- india %>%
  filter(state== "Haryana")
HimachalPradesh <- india %>%
  filter(state== "Himachal Pradesh")
#JammuandKashmir <- india %>%
  #filter(state== "Jammu and Kashmir")
Jharkhand <- india %>%
  filter(state== "Jharkhand")
Karnataka <- india %>%
  filter(state== "Karnataka")
Kerala <- india %>%
  filter(state== "Kerala")
#Lakshadweep <- india %>%
  #filter(state== "Lakshadweep")
MadhyaPradesh <- india %>%
  filter(state== "Madhya Pradesh")
Maharashtra<- india %>%
  filter(state== "Maharashtra")
Manipur <- india %>%
  filter(state== "Manipur")
Meghalaya <- india %>%
  filter(state== "Meghalaya")
Mizoram <- india %>%
  filter(state== "Mizoram")
Nagaland <- india %>%
  filter(state== "Nagaland")
Odisha <- india %>%
  filter(state== "Odisha")
Puducherry <- india %>%
  filter(state== "Puducherry")
Punjab <- india %>%
  filter(state== "Punjab")
Rajasthan <- india %>%
  filter(state== "Rajasthan")
Sikkim <- india %>%
  filter(state== "Sikkim")
TamilNadu <- india %>%
  filter(state== "Tamil Nadu")
Telangana <- india %>%
  filter(state== "Telangana")
Tripura <- india %>%
  filter(state== "Tripura")
UttarPradesh <- india %>%
  filter(state== "Uttar Pradesh")
Uttarakhand<- india %>%
  filter(state== "Uttarakhand")
WestBengal <- india %>%
  filter(state== "West Bengal")



#DIABETES (not separated by rural_urban)
diab_hh_aani <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data= AndamanandNicobarIslands)
calc.icc(diab_hh_aani)
diab_hh_ap <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=AndhraPradesh)
calc.icc(diab_hh_ap)
diab_hh_arp <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=ArunachalPradesh)
calc.icc(diab_hh_arp)
diab_hh_as <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Assam)
calc.icc(diab_hh_as)
diab_hh_bi <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Bihar)
calc.icc(diab_hh_bi)
diab_hh_ch <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Chandigarh)
calc.icc(diab_hh_ch)
diab_hh_chh <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Chhattisgarh)
calc.icc(diab_hh_chh)
#diab_hh_danh <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=DadraandNagarHaveli)
#calc.icc(diab_hh_danh)
diab_hh_dad <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=DamanandDiu)
calc.icc(diab_hh_dad)
diab_hh_de <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Delhi)
calc.icc(diab_hh_de)
diab_hh_goa <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Goa)
calc.icc(diab_hh_goa)
#diab_hh_gu <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Gujarat)
#calc.icc(diab_hh_gu)
diab_hh_ha <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Haryana)
calc.icc(diab_hh_ha)
diab_hh_hp <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=HimachalPradesh)
calc.icc(diab_hh_hp)
#diab_hh_jk <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=JammuandKashmir)
#calc.icc(diab_hh_jk)
diab_hh_jh <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Jharkhand)
calc.icc(diab_hh_jh)
diab_hh_ka <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Karnataka)
calc.icc(diab_hh_ka)
diab_hh_ke <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Kerala)
calc.icc(diab_hh_ke)
#diab_hh_la <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Lakshadweep)
#calc.icc(diab_hh_la)
diab_hh_mp <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=MadhyaPradesh)
calc.icc(diab_hh_mp)
diab_hh_mah <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Maharashtra)
calc.icc(diab_hh_mah)
diab_hh_man <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Manipur)
calc.icc(diab_hh_man)
diab_hh_me <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Meghalaya)
calc.icc(diab_hh_me)
diab_hh_mi <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Mizoram)
calc.icc(diab_hh_mi)
diab_hh_na <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Nagaland)
calc.icc(diab_hh_na)
diab_hh_od <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Odisha)
calc.icc(diab_hh_od)
diab_hh_pud <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Puducherry)
calc.icc(diab_hh_pud)
diab_hh_pun <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Punjab)
calc.icc(diab_hh_pun)
diab_hh_ra <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Rajasthan)
calc.icc(diab_hh_ra)
diab_hh_si <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Sikkim)
calc.icc(diab_hh_si)
diab_hh_tn <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=TamilNadu)
calc.icc(diab_hh_tn)
diab_hh_te <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Telangana)
calc.icc(diab_hh_te)
diab_hh_tri <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Tripura)
calc.icc(diab_hh_tri)
diab_hh_up <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=UttarPradesh)
calc.icc(diab_hh_up)
diab_hh_ut <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=Uttarakhand)
calc.icc(diab_hh_ut)
diab_hh_wb <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=WestBengal)
calc.icc(diab_hh_wb)

htn_hh_aani <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data= AndamanandNicobarIslands)
calc.icc(htn_hh_aani)
htn_hh_ap <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=AndhraPradesh)
calc.icc(htn_hh_ap)
htn_hh_arp <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=ArunachalPradesh)
calc.icc(htn_hh_arp)
htn_hh_as <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Assam)
calc.icc(htn_hh_as)
htn_hh_bi <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Bihar)
calc.icc(htn_hh_bi)
htn_hh_ch <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Chandigarh)
calc.icc(htn_hh_ch)
htn_hh_chh <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Chhattisgarh)
calc.icc(htn_hh_chh)
#htn_hh_danh <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=DadraandNagarHaveli)
#calc.icc(htn_hh_danh)
htn_hh_dad <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=DamanandDiu)
calc.icc(htn_hh_dad)
htn_hh_de <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Delhi)
calc.icc(htn_hh_de)
htn_hh_goa <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Goa)
calc.icc(htn_hh_goa)
#htn_hh_gu <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Gujarat)
#calc.icc(htn_hh_gu)
htn_hh_ha <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Haryana)
calc.icc(htn_hh_ha)
htn_hh_hp <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=HimachalPradesh)
calc.icc(htn_hh_hp)
#htn_hh_jk <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=JammuandKashmir)
#calc.icc(htn_hh_jk)
htn_hh_jh <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Jharkhand)
calc.icc(htn_hh_jh)
htn_hh_ka <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Karnataka)
calc.icc(htn_hh_ka)
htn_hh_ke <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Kerala)
calc.icc(htn_hh_ke)
#htn_hh_la <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Lakshadweep)
#calc.icc(htn_hh_la)
htn_hh_mp <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=MadhyaPradesh)
calc.icc(htn_hh_mp)
htn_hh_mah <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Maharashtra)
calc.icc(htn_hh_mah)
htn_hh_man <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Manipur)
calc.icc(htn_hh_man)
htn_hh_me <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Meghalaya)
calc.icc(htn_hh_me)
htn_hh_mi <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Mizoram)
calc.icc(htn_hh_mi)
htn_hh_na <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Nagaland)
calc.icc(htn_hh_na)
htn_hh_od <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Odisha)
calc.icc(htn_hh_od)
htn_hh_pud <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Puducherry)
calc.icc(htn_hh_pud)
htn_hh_pun <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Punjab)
calc.icc(htn_hh_pun)
htn_hh_ra <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Rajasthan)
calc.icc(htn_hh_ra)
htn_hh_si <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Sikkim)
calc.icc(htn_hh_si)
htn_hh_tn <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=TamilNadu)
calc.icc(htn_hh_tn)
htn_hh_te <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Telangana)
calc.icc(htn_hh_te)
htn_hh_tri <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Tripura)
calc.icc(htn_hh_tri)
htn_hh_up <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=UttarPradesh)
calc.icc(htn_hh_up)
htn_hh_ut <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=Uttarakhand)
calc.icc(htn_hh_ut)
htn_hh_wb <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=WestBengal)
calc.icc(htn_hh_wb)


csmoke_hh_aani <- lmer(csmoke ~ 1 + (1|hh_id), data= AndamanandNicobarIslands)
calc.icc(csmoke_hh_aani)
csmoke_hh_ap <- lmer(csmoke ~ 1 + (1|hh_id), data=AndhraPradesh)
calc.icc(csmoke_hh_ap)
csmoke_hh_arp <- lmer(csmoke ~ 1 + (1|hh_id), data=ArunachalPradesh)
calc.icc(csmoke_hh_arp)
csmoke_hh_as <- lmer(csmoke ~ 1 + (1|hh_id), data=Assam)
calc.icc(csmoke_hh_as)
csmoke_hh_bi <- lmer(csmoke ~ 1 + (1|hh_id), data=Bihar)
calc.icc(csmoke_hh_bi)
csmoke_hh_ch <- lmer(csmoke ~ 1 + (1|hh_id), data=Chandigarh)
calc.icc(csmoke_hh_ch)
csmoke_hh_chh <- lmer(csmoke ~ 1 + (1|hh_id), data=Chhattisgarh)
calc.icc(csmoke_hh_chh)
#csmoke_hh_danh <- lmer(csmoke ~ 1 + (1|hh_id), data=DadraandNagarHaveli)
#calc.icc(csmoke_hh_danh)
csmoke_hh_dad <- lmer(csmoke ~ 1 + (1|hh_id), data=DamanandDiu)
calc.icc(csmoke_hh_dad)
csmoke_hh_de <- lmer(csmoke ~ 1 + (1|hh_id), data=Delhi)
calc.icc(csmoke_hh_de)
csmoke_hh_goa <- lmer(csmoke ~ 1 + (1|hh_id), data=Goa)
calc.icc(csmoke_hh_goa)
#csmoke_hh_gu <- lmer(csmoke ~ 1 + (1|hh_id), data=Gujarat)
#calc.icc(csmoke_hh_gu)
csmoke_hh_ha <- lmer(csmoke ~ 1 + (1|hh_id), data=Haryana)
calc.icc(csmoke_hh_ha)
csmoke_hh_hp <- lmer(csmoke ~ 1 + (1|hh_id), data=HimachalPradesh)
calc.icc(csmoke_hh_hp)
#csmoke_hh_jk <- lmer(csmoke ~ 1 + (1|hh_id), data=JammuandKashmir)
#calc.icc(csmoke_hh_jk)
csmoke_hh_jh <- lmer(csmoke ~ 1 + (1|hh_id), data=Jharkhand)
calc.icc(csmoke_hh_jh)
csmoke_hh_ka <- lmer(csmoke ~ 1 + (1|hh_id), data=Karnataka)
calc.icc(csmoke_hh_ka)
csmoke_hh_ke <- lmer(csmoke ~ 1 + (1|hh_id), data=Kerala)
calc.icc(csmoke_hh_ke)
#csmoke_hh_la <- lmer(csmoke ~ 1 + (1|hh_id), data=Lakshadweep)
#calc.icc(csmoke_hh_la)
csmoke_hh_mp <- lmer(csmoke ~ 1 + (1|hh_id), data=MadhyaPradesh)
calc.icc(csmoke_hh_mp)
csmoke_hh_mah <- lmer(csmoke ~ 1 + (1|hh_id), data=Maharashtra)
calc.icc(csmoke_hh_mah)
csmoke_hh_man <- lmer(csmoke ~ 1 + (1|hh_id), data=Manipur)
calc.icc(csmoke_hh_man)
csmoke_hh_me <- lmer(csmoke ~ 1 + (1|hh_id), data=Meghalaya)
calc.icc(csmoke_hh_me)
csmoke_hh_mi <- lmer(csmoke ~ 1 + (1|hh_id), data=Mizoram)
calc.icc(csmoke_hh_mi)
csmoke_hh_na <- lmer(csmoke ~ 1 + (1|hh_id), data=Nagaland)
calc.icc(csmoke_hh_na)
csmoke_hh_od <- lmer(csmoke ~ 1 + (1|hh_id), data=Odisha)
calc.icc(csmoke_hh_od)
csmoke_hh_pud <- lmer(csmoke ~ 1 + (1|hh_id), data=Puducherry)
calc.icc(csmoke_hh_pud)
csmoke_hh_pun <- lmer(csmoke ~ 1 + (1|hh_id), data=Punjab)
calc.icc(csmoke_hh_pun)
csmoke_hh_ra <- lmer(csmoke ~ 1 + (1|hh_id), data=Rajasthan)
calc.icc(csmoke_hh_ra)
csmoke_hh_si <- lmer(csmoke ~ 1 + (1|hh_id), data=Sikkim)
calc.icc(csmoke_hh_si)
csmoke_hh_tn <- lmer(csmoke ~ 1 + (1|hh_id), data=TamilNadu)
calc.icc(csmoke_hh_tn)
csmoke_hh_te <- lmer(csmoke ~ 1 + (1|hh_id), data=Telangana)
calc.icc(csmoke_hh_te)
csmoke_hh_tri <- lmer(csmoke ~ 1 + (1|hh_id), data=Tripura)
calc.icc(csmoke_hh_tri)
csmoke_hh_up <- lmer(csmoke ~ 1 + (1|hh_id), data=UttarPradesh)
calc.icc(csmoke_hh_up)
csmoke_hh_ut <- lmer(csmoke ~ 1 + (1|hh_id), data=Uttarakhand)
calc.icc(csmoke_hh_ut)
csmoke_hh_wb <- lmer(csmoke ~ 1 + (1|hh_id), data=WestBengal)
calc.icc(csmoke_hh_wb)

overweight_hh_aani <- lmer(overweight23 ~ 1 + (1|hh_id), data= AndamanandNicobarIslands)
calc.icc(overweight_hh_aani)
overweight_hh_ap <- lmer(overweight23 ~ 1 + (1|hh_id), data=AndhraPradesh)
calc.icc(overweight_hh_ap)
overweight_hh_arp <- lmer(overweight23 ~ 1 + (1|hh_id), data=ArunachalPradesh)
calc.icc(overweight_hh_arp)
overweight_hh_as <- lmer(overweight23 ~ 1 + (1|hh_id), data=Assam)
calc.icc(overweight_hh_as)
overweight_hh_bi <- lmer(overweight23 ~ 1 + (1|hh_id), data=Bihar)
calc.icc(overweight_hh_bi)
overweight_hh_ch <- lmer(overweight23 ~ 1 + (1|hh_id), data=Chandigarh)
calc.icc(overweight_hh_ch)
overweight_hh_chh <- lmer(overweight23 ~ 1 + (1|hh_id), data=Chhattisgarh)
calc.icc(overweight_hh_chh)
#overweight_hh_danh <- lmer(overweight23 ~ 1 + (1|hh_id), data=DadraandNagarHaveli)
#calc.icc(overweight_hh_danh)
overweight_hh_dad <- lmer(overweight23 ~ 1 + (1|hh_id), data=DamanandDiu)
calc.icc(overweight_hh_dad)
overweight_hh_de <- lmer(overweight23 ~ 1 + (1|hh_id), data=Delhi)
calc.icc(overweight_hh_de)
overweight_hh_goa <- lmer(overweight23 ~ 1 + (1|hh_id), data=Goa)
calc.icc(overweight_hh_goa)
#overweight_hh_gu <- lmer(overweight23 ~ 1 + (1|hh_id), data=Gujarat)
#calc.icc(overweight_hh_gu)
overweight_hh_ha <- lmer(overweight23 ~ 1 + (1|hh_id), data=Haryana)
calc.icc(overweight_hh_ha)
overweight_hh_hp <- lmer(overweight23 ~ 1 + (1|hh_id), data=HimachalPradesh)
calc.icc(overweight_hh_hp)
#overweight_hh_jk <- lmer(overweight23 ~ 1 + (1|hh_id), data=JammuandKashmir)
#calc.icc(overweight_hh_jk)
overweight_hh_jh <- lmer(overweight23 ~ 1 + (1|hh_id), data=Jharkhand)
calc.icc(overweight_hh_jh)
overweight_hh_ka <- lmer(overweight23 ~ 1 + (1|hh_id), data=Karnataka)
calc.icc(overweight_hh_ka)
overweight_hh_ke <- lmer(overweight23 ~ 1 + (1|hh_id), data=Kerala)
calc.icc(overweight_hh_ke)
#overweight_hh_la <- lmer(overweight23 ~ 1 + (1|hh_id), data=Lakshadweep)
#calc.icc(overweight_hh_la)
overweight_hh_mp <- lmer(overweight23 ~ 1 + (1|hh_id), data=MadhyaPradesh)
calc.icc(overweight_hh_mp)
overweight_hh_mah <- lmer(overweight23 ~ 1 + (1|hh_id), data=Maharashtra)
calc.icc(overweight_hh_mah)
overweight_hh_man <- lmer(overweight23 ~ 1 + (1|hh_id), data=Manipur)
calc.icc(overweight_hh_man)
overweight_hh_me <- lmer(overweight23 ~ 1 + (1|hh_id), data=Meghalaya)
calc.icc(overweight_hh_me)
overweight_hh_mi <- lmer(overweight23 ~ 1 + (1|hh_id), data=Mizoram)
calc.icc(overweight_hh_mi)
overweight_hh_na <- lmer(overweight23 ~ 1 + (1|hh_id), data=Nagaland)
calc.icc(overweight_hh_na)
overweight_hh_od <- lmer(overweight23 ~ 1 + (1|hh_id), data=Odisha)
calc.icc(overweight_hh_od)
overweight_hh_pud <- lmer(overweight23 ~ 1 + (1|hh_id), data=Puducherry)
calc.icc(overweight_hh_pud)
overweight_hh_pun <- lmer(overweight23 ~ 1 + (1|hh_id), data=Punjab)
calc.icc(overweight_hh_pun)
overweight_hh_ra <- lmer(overweight23 ~ 1 + (1|hh_id), data=Rajasthan)
calc.icc(overweight_hh_ra)
overweight_hh_si <- lmer(overweight23 ~ 1 + (1|hh_id), data=Sikkim)
calc.icc(overweight_hh_si)
overweight_hh_tn <- lmer(overweight23 ~ 1 + (1|hh_id), data=TamilNadu)
calc.icc(overweight_hh_tn)
overweight_hh_te <- lmer(overweight23 ~ 1 + (1|hh_id), data=Telangana)
calc.icc(overweight_hh_te)
overweight_hh_tri <- lmer(overweight23 ~ 1 + (1|hh_id), data=Tripura)
calc.icc(overweight_hh_tri)
overweight_hh_up <- lmer(overweight23 ~ 1 + (1|hh_id), data=UttarPradesh)
calc.icc(overweight_hh_up)
overweight_hh_ut <- lmer(overweight23 ~ 1 + (1|hh_id), data=Uttarakhand)
calc.icc(overweight_hh_ut)
overweight_hh_wb <- lmer(overweight23 ~ 1 + (1|hh_id), data=WestBengal)
calc.icc(overweight_hh_wb)


obese27_hh_aani <- lmer(obese27 ~ 1 + (1|hh_id), data= AndamanandNicobarIslands)
calc.icc(obese27_hh_aani)
obese27_hh_ap <- lmer(obese27 ~ 1 + (1|hh_id), data=AndhraPradesh)
calc.icc(obese27_hh_ap)
obese27_hh_arp <- lmer(obese27 ~ 1 + (1|hh_id), data=ArunachalPradesh)
calc.icc(obese27_hh_arp)
obese27_hh_as <- lmer(obese27 ~ 1 + (1|hh_id), data=Assam)
calc.icc(obese27_hh_as)
obese27_hh_bi <- lmer(obese27 ~ 1 + (1|hh_id), data=Bihar)
calc.icc(obese27_hh_bi)
obese27_hh_ch <- lmer(obese27 ~ 1 + (1|hh_id), data=Chandigarh)
calc.icc(obese27_hh_ch)
obese27_hh_chh <- lmer(obese27 ~ 1 + (1|hh_id), data=Chhattisgarh)
calc.icc(obese27_hh_chh)
#obese27_hh_danh <- lmer(obese27 ~ 1 + (1|hh_id), data=DadraandNagarHaveli)
#calc.icc(obese27_hh_danh)
obese27_hh_dad <- lmer(obese27 ~ 1 + (1|hh_id), data=DamanandDiu)
calc.icc(obese27_hh_dad)
obese27_hh_de <- lmer(obese27 ~ 1 + (1|hh_id), data=Delhi)
calc.icc(obese27_hh_de)
obese27_hh_goa <- lmer(obese27 ~ 1 + (1|hh_id), data=Goa)
calc.icc(obese27_hh_goa)
#obese27_hh_gu <- lmer(obese27 ~ 1 + (1|hh_id), data=Gujarat)
#calc.icc(obese27_hh_gu)
obese27_hh_ha <- lmer(obese27 ~ 1 + (1|hh_id), data=Haryana)
calc.icc(obese27_hh_ha)
obese27_hh_hp <- lmer(obese27 ~ 1 + (1|hh_id), data=HimachalPradesh)
calc.icc(obese27_hh_hp)
#obese27_hh_jk <- lmer(obese27 ~ 1 + (1|hh_id), data=JammuandKashmir)
#calc.icc(obese27_hh_jk)
obese27_hh_jh <- lmer(obese27 ~ 1 + (1|hh_id), data=Jharkhand)
calc.icc(obese27_hh_jh)
obese27_hh_ka <- lmer(obese27 ~ 1 + (1|hh_id), data=Karnataka)
calc.icc(obese27_hh_ka)
obese27_hh_ke <- lmer(obese27 ~ 1 + (1|hh_id), data=Kerala)
calc.icc(obese27_hh_ke)
#obese27_hh_la <- lmer(obese27 ~ 1 + (1|hh_id), data=Lakshadweep)
#calc.icc(obese27_hh_la)
obese27_hh_mp <- lmer(obese27 ~ 1 + (1|hh_id), data=MadhyaPradesh)
calc.icc(obese27_hh_mp)
obese27_hh_mah <- lmer(obese27 ~ 1 + (1|hh_id), data=Maharashtra)
calc.icc(obese27_hh_mah)
obese27_hh_man <- lmer(obese27 ~ 1 + (1|hh_id), data=Manipur)
calc.icc(obese27_hh_man)
obese27_hh_me <- lmer(obese27 ~ 1 + (1|hh_id), data=Meghalaya)
calc.icc(obese27_hh_me)
obese27_hh_mi <- lmer(obese27 ~ 1 + (1|hh_id), data=Mizoram)
calc.icc(obese27_hh_mi)
obese27_hh_na <- lmer(obese27 ~ 1 + (1|hh_id), data=Nagaland)
calc.icc(obese27_hh_na)
obese27_hh_od <- lmer(obese27 ~ 1 + (1|hh_id), data=Odisha)
calc.icc(obese27_hh_od)
obese27_hh_pud <- lmer(obese27 ~ 1 + (1|hh_id), data=Puducherry)
calc.icc(obese27_hh_pud)
obese27_hh_pun <- lmer(obese27 ~ 1 + (1|hh_id), data=Punjab)
calc.icc(obese27_hh_pun)
obese27_hh_ra <- lmer(obese27 ~ 1 + (1|hh_id), data=Rajasthan)
calc.icc(obese27_hh_ra)
obese27_hh_si <- lmer(obese27 ~ 1 + (1|hh_id), data=Sikkim)
calc.icc(obese27_hh_si)
obese27_hh_tn <- lmer(obese27 ~ 1 + (1|hh_id), data=TamilNadu)
calc.icc(obese27_hh_tn)
obese27_hh_te <- lmer(obese27 ~ 1 + (1|hh_id), data=Telangana)
calc.icc(obese27_hh_te)
obese27_hh_tri <- lmer(obese27 ~ 1 + (1|hh_id), data=Tripura)
calc.icc(obese27_hh_tri)
obese27_hh_up <- lmer(obese27 ~ 1 + (1|hh_id), data=UttarPradesh)
calc.icc(obese27_hh_up)
obese27_hh_ut <- lmer(obese27 ~ 1 + (1|hh_id), data=Uttarakhand)
calc.icc(obese27_hh_ut)
obese27_hh_wb <- lmer(obese27 ~ 1 + (1|hh_id), data=WestBengal)
calc.icc(obese27_hh_wb)



############calculate bootstraps for state iccs########

boot_diab_hh_aani <- bootMer(diab_hh_aani, calc.icc, nsim=500)
boot_diab_hh_ap <- bootMer(diab_hh_ap, calc.icc, nsim=500)
boot_diab_hh_arp <- bootMer(diab_hh_arp, calc.icc, nsim=500)
boot_diab_hh_as <- bootMer(diab_hh_as, calc.icc, nsim=500)
boot_diab_hh_bi <- bootMer(diab_hh_bi, calc.icc, nsim=500)
boot_diab_hh_ch <- bootMer(diab_hh_ch, calc.icc, nsim=500)
boot_diab_hh_chh <- bootMer(diab_hh_chh, calc.icc, nsim=500)
#boot_diab_hh_danh <- bootMer(diab_hh_danh, calc.icc, nsim=500)
boot_diab_hh_dad <- bootMer(diab_hh_dad, calc.icc, nsim=500)
boot_diab_hh_de <- bootMer(diab_hh_de, calc.icc, nsim=500) 
boot_diab_hh_goa <- bootMer(diab_hh_goa, calc.icc, nsim=500)
#boot_diab_hh_gu <- bootMer(diab_hh_gu, calc.icc, nsim=500)
boot_diab_hh_ha <- bootMer(diab_hh_ha, calc.icc, nsim=500)
boot_diab_hh_hp <- bootMer(diab_hh_hp, calc.icc, nsim=500)
#boot_diab_hh_jk <- bootMer(diab_hh_jk, calc.icc, nsim=500)
boot_diab_hh_jh <- bootMer(diab_hh_jh, calc.icc, nsim=500)
boot_diab_hh_ka <- bootMer(diab_hh_ka, calc.icc, nsim=500)
boot_diab_hh_ke <- bootMer(diab_hh_ke, calc.icc, nsim=500)
#boot_diab_hh_la <- bootMer(diab_hh_la, calc.icc, nsim=500)
boot_diab_hh_mp <- bootMer(diab_hh_mp, calc.icc, nsim=500)
boot_diab_hh_mah <- bootMer(diab_hh_mah, calc.icc, nsim=500)
boot_diab_hh_man <- bootMer(diab_hh_man, calc.icc, nsim=500)
boot_diab_hh_me <- bootMer(diab_hh_me, calc.icc, nsim=500)
boot_diab_hh_mi <- bootMer(diab_hh_mi, calc.icc, nsim=500) 
boot_diab_hh_na <- bootMer(diab_hh_na, calc.icc, nsim=500)
boot_diab_hh_od <- bootMer(diab_hh_od, calc.icc, nsim=500)
boot_diab_hh_pud <- bootMer(diab_hh_pud, calc.icc, nsim=500)
boot_diab_hh_pun <- bootMer(diab_hh_pun, calc.icc, nsim=500)
boot_diab_hh_ra <- bootMer(diab_hh_ra, calc.icc, nsim=500)
boot_diab_hh_si <- bootMer(diab_hh_si, calc.icc, nsim=500)
boot_diab_hh_tn <- bootMer(diab_hh_tn, calc.icc, nsim=500)
boot_diab_hh_te <- bootMer(diab_hh_te, calc.icc, nsim=500)
boot_diab_hh_tri <- bootMer(diab_hh_tri, calc.icc, nsim=500)
boot_diab_hh_up <- bootMer(diab_hh_up, calc.icc, nsim=500)
boot_diab_hh_ut <- bootMer(diab_hh_ut, calc.icc, nsim=500)
boot_diab_hh_wb <- bootMer(diab_hh_wb, calc.icc, nsim=500)

boot_htn_hh_aani <- bootMer(htn_hh_aani, calc.icc, nsim=500)
boot_htn_hh_ap <- bootMer(htn_hh_ap, calc.icc, nsim=500)
boot_htn_hh_arp <- bootMer(htn_hh_arp, calc.icc, nsim=500)
boot_htn_hh_as <- bootMer(htn_hh_as, calc.icc, nsim=500)
boot_htn_hh_bi <- bootMer(htn_hh_bi, calc.icc, nsim=500)
boot_htn_hh_ch <- bootMer(htn_hh_ch, calc.icc, nsim=500)
boot_htn_hh_chh <- bootMer(htn_hh_chh, calc.icc, nsim=500)
#boot_htn_hh_danh <- bootMer(htn_hh_danh, calc.icc, nsim=500)
boot_htn_hh_dad <- bootMer(htn_hh_dad, calc.icc, nsim=500)
boot_htn_hh_de <- bootMer(htn_hh_de, calc.icc, nsim=500)
boot_htn_hh_goa <- bootMer(htn_hh_goa, calc.icc, nsim=500)
#boot_htn_hh_gu <- bootMer(htn_hh_gu, calc.icc, nsim=500)
boot_htn_hh_ha <- bootMer(htn_hh_ha, calc.icc, nsim=500)
boot_htn_hh_hp <- bootMer(htn_hh_hp, calc.icc, nsim=500)
#boot_htn_hh_jk <- bootMer(htn_hh_jk, calc.icc, nsim=500)
boot_htn_hh_jh <- bootMer(htn_hh_jh, calc.icc, nsim=500)
boot_htn_hh_ka <- bootMer(htn_hh_ka, calc.icc, nsim=500)
boot_htn_hh_ke <- bootMer(htn_hh_ke, calc.icc, nsim=500)
#boot_htn_hh_la <- bootMer(htn_hh_la, calc.icc, nsim=500)
boot_htn_hh_mp <- bootMer(htn_hh_mp, calc.icc, nsim=500)
boot_htn_hh_mah <- bootMer(htn_hh_mah, calc.icc, nsim=500)
boot_htn_hh_man <- bootMer(htn_hh_man, calc.icc, nsim=500)
boot_htn_hh_me <- bootMer(htn_hh_me, calc.icc, nsim=500)
boot_htn_hh_mi <- bootMer(htn_hh_mi, calc.icc, nsim=500)
boot_htn_hh_na <- bootMer(htn_hh_na, calc.icc, nsim=500)
boot_htn_hh_od <- bootMer(htn_hh_od, calc.icc, nsim=500)
boot_htn_hh_pud <- bootMer(htn_hh_pud, calc.icc, nsim=500)
boot_htn_hh_pun <- bootMer(htn_hh_pun, calc.icc, nsim=500)
boot_htn_hh_ra <- bootMer(htn_hh_ra, calc.icc, nsim=500)
boot_htn_hh_si <- bootMer(htn_hh_si, calc.icc, nsim=500)
boot_htn_hh_tn <- bootMer(htn_hh_tn, calc.icc, nsim=500)
boot_htn_hh_te <- bootMer(htn_hh_te, calc.icc, nsim=500)
boot_htn_hh_tri <- bootMer(htn_hh_tri, calc.icc, nsim=500)
boot_htn_hh_up <- bootMer(htn_hh_up, calc.icc, nsim=500)
boot_htn_hh_ut <- bootMer(htn_hh_ut, calc.icc, nsim=500)
boot_htn_hh_wb <- bootMer(htn_hh_wb, calc.icc, nsim=500)

boot_csmoke_hh_aani <- bootMer(csmoke_hh_aani, calc.icc, nsim=500)
boot_csmoke_hh_ap <- bootMer(csmoke_hh_ap, calc.icc, nsim=500)
boot_csmoke_hh_arp <- bootMer(csmoke_hh_arp, calc.icc, nsim=500)
boot_csmoke_hh_as <- bootMer(csmoke_hh_as, calc.icc, nsim=500)
boot_csmoke_hh_bi <- bootMer(csmoke_hh_bi, calc.icc, nsim=500)
boot_csmoke_hh_ch <- bootMer(csmoke_hh_ch, calc.icc, nsim=500)
boot_csmoke_hh_chh <- bootMer(csmoke_hh_chh, calc.icc, nsim=500)
#boot_csmoke_hh_danh <- bootMer(csmoke_hh_danh, calc.icc, nsim=500)
boot_csmoke_hh_dad <- bootMer(csmoke_hh_dad, calc.icc, nsim=500)
boot_csmoke_hh_de <- bootMer(csmoke_hh_de, calc.icc, nsim=500)
boot_csmoke_hh_goa <- bootMer(csmoke_hh_goa, calc.icc, nsim=500)
#boot_csmoke_hh_gu <- bootMer(csmoke_hh_gu, calc.icc, nsim=500)
boot_csmoke_hh_ha <- bootMer(csmoke_hh_ha, calc.icc, nsim=500)
boot_csmoke_hh_hp <- bootMer(csmoke_hh_hp, calc.icc, nsim=500)
#boot_csmoke_hh_jk <- bootMer(csmoke_hh_jk, calc.icc, nsim=500)
boot_csmoke_hh_jh <- bootMer(csmoke_hh_jh, calc.icc, nsim=500)
boot_csmoke_hh_ka <- bootMer(csmoke_hh_ka, calc.icc, nsim=500)
boot_csmoke_hh_ke <- bootMer(csmoke_hh_ke, calc.icc, nsim=500)
#boot_csmoke_hh_la <- bootMer(csmoke_hh_la, calc.icc, nsim=500)
boot_csmoke_hh_mp <- bootMer(csmoke_hh_mp, calc.icc, nsim=500)
boot_csmoke_hh_mah <- bootMer(csmoke_hh_mah, calc.icc, nsim=500)
boot_csmoke_hh_man <- bootMer(csmoke_hh_man, calc.icc, nsim=500)
boot_csmoke_hh_me <- bootMer(csmoke_hh_me, calc.icc, nsim=500)
boot_csmoke_hh_mi <- bootMer(csmoke_hh_mi, calc.icc, nsim=500)
boot_csmoke_hh_na <- bootMer(csmoke_hh_na, calc.icc, nsim=500)
boot_csmoke_hh_od <- bootMer(csmoke_hh_od, calc.icc, nsim=500)
boot_csmoke_hh_pud <- bootMer(csmoke_hh_pud, calc.icc, nsim=500)
boot_csmoke_hh_pun <- bootMer(csmoke_hh_pun, calc.icc, nsim=500)
boot_csmoke_hh_ra <- bootMer(csmoke_hh_ra, calc.icc, nsim=500)
boot_csmoke_hh_si <- bootMer(csmoke_hh_si, calc.icc, nsim=500)
boot_csmoke_hh_tn <- bootMer(csmoke_hh_tn, calc.icc, nsim=500)
boot_csmoke_hh_te <- bootMer(csmoke_hh_te, calc.icc, nsim=500)
boot_csmoke_hh_tri <- bootMer(csmoke_hh_tri, calc.icc, nsim=500)
boot_csmoke_hh_up <- bootMer(csmoke_hh_up, calc.icc, nsim=500)
boot_csmoke_hh_ut <- bootMer(csmoke_hh_ut, calc.icc, nsim=500)
boot_csmoke_hh_wb <- bootMer(csmoke_hh_wb, calc.icc, nsim=500)

boot_overweight_hh_aani <- bootMer(overweight_hh_aani, calc.icc, nsim=500)
boot_overweight_hh_ap <- bootMer(overweight_hh_ap, calc.icc, nsim=500)
boot_overweight_hh_arp <- bootMer(overweight_hh_arp, calc.icc, nsim=500)
boot_overweight_hh_as <- bootMer(overweight_hh_as, calc.icc, nsim=500)
boot_overweight_hh_bi <- bootMer(overweight_hh_bi, calc.icc, nsim=500)
boot_overweight_hh_ch <- bootMer(overweight_hh_ch, calc.icc, nsim=500)
boot_overweight_hh_chh <- bootMer(overweight_hh_chh, calc.icc, nsim=500)
#boot_overweight_hh_danh <- bootMer(overweight_hh_danh, calc.icc, nsim=500)
boot_overweight_hh_dad <- bootMer(overweight_hh_dad, calc.icc, nsim=500)
boot_overweight_hh_de <- bootMer(overweight_hh_de, calc.icc, nsim=500)
boot_overweight_hh_goa <- bootMer(overweight_hh_goa, calc.icc, nsim=500)
#boot_overweight_hh_gu <- bootMer(overweight_hh_gu, calc.icc, nsim=500)
boot_overweight_hh_ha <- bootMer(overweight_hh_ha, calc.icc, nsim=500)
boot_overweight_hh_hp <- bootMer(overweight_hh_hp, calc.icc, nsim=500)
#boot_overweight_hh_jk <- bootMer(overweight_hh_jk, calc.icc, nsim=500)
boot_overweight_hh_jh <- bootMer(overweight_hh_jh, calc.icc, nsim=500)
boot_overweight_hh_ka <- bootMer(overweight_hh_ka, calc.icc, nsim=500)
boot_overweight_hh_ke <- bootMer(overweight_hh_ke, calc.icc, nsim=500)
#boot_overweight_hh_la <- bootMer(overweight_hh_la, calc.icc, nsim=500)
boot_overweight_hh_mp <- bootMer(overweight_hh_mp, calc.icc, nsim=500)
boot_overweight_hh_mah <- bootMer(overweight_hh_mah, calc.icc, nsim=500)
boot_overweight_hh_man <- bootMer(overweight_hh_man, calc.icc, nsim=500)
boot_overweight_hh_me <- bootMer(overweight_hh_me, calc.icc, nsim=500)
boot_overweight_hh_mi <- bootMer(overweight_hh_mi, calc.icc, nsim=500)

boot_overweight_hh_na <- bootMer(overweight_hh_na, calc.icc, nsim=500)
boot_overweight_hh_od <- bootMer(overweight_hh_od, calc.icc, nsim=500)
boot_overweight_hh_pud <- bootMer(overweight_hh_pud, calc.icc, nsim=500)
boot_overweight_hh_pun <- bootMer(overweight_hh_pun, calc.icc, nsim=500)
boot_overweight_hh_ra <- bootMer(overweight_hh_ra, calc.icc, nsim=500)
boot_overweight_hh_si <- bootMer(overweight_hh_si, calc.icc, nsim=500)
boot_overweight_hh_tn <- bootMer(overweight_hh_tn, calc.icc, nsim=500)
boot_overweight_hh_te <- bootMer(overweight_hh_te, calc.icc, nsim=500)
boot_overweight_hh_tri <- bootMer(overweight_hh_tri, calc.icc, nsim=500)
boot_overweight_hh_up <- bootMer(overweight_hh_up, calc.icc, nsim=500)
boot_overweight_hh_ut <- bootMer(overweight_hh_ut, calc.icc, nsim=500)
boot_overweight_hh_wb <- bootMer(overweight_hh_wb, calc.icc, nsim=500)

boot_obese27_hh_aani <- bootMer(obese27_hh_aani, calc.icc, nsim=500)
boot_obese27_hh_ap <- bootMer(obese27_hh_ap, calc.icc, nsim=500)
boot_obese27_hh_arp <- bootMer(obese27_hh_arp, calc.icc, nsim=500)
boot_obese27_hh_as <- bootMer(obese27_hh_as, calc.icc, nsim=500)
boot_obese27_hh_bi <- bootMer(obese27_hh_bi, calc.icc, nsim=500)
boot_obese27_hh_ch <- bootMer(obese27_hh_ch, calc.icc, nsim=500)
boot_obese27_hh_chh <- bootMer(obese27_hh_chh, calc.icc, nsim=500)
#boot_obese27_hh_danh <- bootMer(obese27_hh_danh, calc.icc, nsim=500)
boot_obese27_hh_dad <- bootMer(obese27_hh_dad, calc.icc, nsim=500)
boot_obese27_hh_de <- bootMer(obese27_hh_de, calc.icc, nsim=500)
boot_obese27_hh_goa <- bootMer(obese27_hh_goa, calc.icc, nsim=500)
#boot_obese27_hh_gu <- bootMer(obese27_hh_gu, calc.icc, nsim=500)
boot_obese27_hh_ha <- bootMer(obese27_hh_ha, calc.icc, nsim=500)
boot_obese27_hh_hp <- bootMer(obese27_hh_hp, calc.icc, nsim=500)
#boot_obese27_hh_jk <- bootMer(obese27_hh_jk ,calc.icc, nsim=500)
boot_obese27_hh_jh <- bootMer(obese27_hh_jh, calc.icc, nsim=500)
boot_obese27_hh_ka <- bootMer(obese27_hh_ka, calc.icc, nsim=500)
boot_obese27_hh_ke <- bootMer(obese27_hh_ke, calc.icc, nsim=500)
#boot_obese27_hh_la <- bootMer(obese27_hh_la, calc.icc, nsim=500)
boot_obese27_hh_mp <- bootMer(obese27_hh_mp, calc.icc, nsim=500)
boot_obese27_hh_mah <- bootMer(obese27_hh_mah, calc.icc, nsim=500)
boot_obese27_hh_man <- bootMer(obese27_hh_man, calc.icc, nsim=500)
boot_obese27_hh_me <- bootMer(obese27_hh_me, calc.icc, nsim=500)
boot_obese27_hh_mi <- bootMer(obese27_hh_mi, calc.icc, nsim=500)
boot_obese27_hh_na <- bootMer(obese27_hh_na, calc.icc, nsim=500)
boot_obese27_hh_od <- bootMer(obese27_hh_od, calc.icc, nsim=500)
boot_obese27_hh_pud <- bootMer(obese27_hh_pud, calc.icc, nsim=500)
boot_obese27_hh_pun <- bootMer(obese27_hh_pun, calc.icc, nsim=500)
boot_obese27_hh_ra <- bootMer(obese27_hh_ra, calc.icc, nsim=500)
boot_obese27_hh_si <- bootMer(obese27_hh_si, calc.icc, nsim=500)
boot_obese27_hh_tn <- bootMer(obese27_hh_tn, calc.icc, nsim=500)
boot_obese27_hh_te <- bootMer(obese27_hh_te, calc.icc, nsim=500)
boot_obese27_hh_tri <- bootMer(obese27_hh_tri, calc.icc, nsim=500)
boot_obese27_hh_up <- bootMer(obese27_hh_up, calc.icc, nsim=500)
boot_obese27_hh_ut <- bootMer(obese27_hh_ut, calc.icc, nsim=500)
boot_obese27_hh_wb <- bootMer(obese27_hh_wb, calc.icc, nsim=500)

############################################

quantile(boot_diab_hh_aani$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_ap$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_arp$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_as$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_bi$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_ch$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_chh$t, c(0.025, 0.975)) 
#quantile(boot_diab_hh_danh$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_dad$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_de$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_goa$t, c(0.025, 0.975)) 
#quantile(boot_diab_hh_gu$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_ha$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_hp$t, c(0.025, 0.975)) 
#quantile(boot_diab_hh_jk$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_jh$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_ka$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_ke$t, c(0.025, 0.975)) 
#quantile(boot_diab_hh_la$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_mp$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_mah$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_man$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_me$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_mi$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_na$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_od$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_pud$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_pun$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_ra$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_si$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_tn$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_te$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_tri$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_up$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_ut$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_wb$t, c(0.025, 0.975)) 

quantile(boot_htn_hh_aani$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_ap$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_arp$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_as$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_bi$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_ch$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_chh$t, c(0.025, 0.975)) 
#quantile(boot_htn_hh_danh$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_dad$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_de$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_goa$t, c(0.025, 0.975)) 
#quantile(boot_htn_hh_gu$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_ha$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_hp$t, c(0.025, 0.975)) 
#quantile(boot_htn_hh_jk$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_jh$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_ka$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_ke$t, c(0.025, 0.975)) 
#quantile(boot_htn_hh_la$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_mp$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_mah$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_man$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_me$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_mi$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_na$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_od$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_pud$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_pun$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_ra$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_si$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_tn$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_te$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_tri$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_up$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_ut$t, c(0.025, 0.975)) 
quantile(boot_htn_hh_wb$t, c(0.025, 0.975)) 

quantile(boot_csmoke_hh_aani$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_ap$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_arp$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_as$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_bi$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_ch$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_chh$t, c(0.025, 0.975)) 
#quantile(boot_csmoke_hh_danh$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_dad$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_de$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_goa$t, c(0.025, 0.975)) 
#quantile(boot_csmoke_hh_gu$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_ha$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_hp$t, c(0.025, 0.975)) 
#quantile(boot_csmoke_hh_jk$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_jh$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_ka$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_ke$t, c(0.025, 0.975)) 
#quantile(boot_csmoke_hh_la$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_mp$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_mah$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_man$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_me$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_mi$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_na$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_od$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_pud$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_pun$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_ra$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_si$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_tn$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_te$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_tri$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_up$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_ut$t, c(0.025, 0.975)) 
quantile(boot_csmoke_hh_wb$t, c(0.025, 0.975)) 

quantile(boot_overweight_hh_aani$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_ap$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_arp$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_as$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_bi$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_ch$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_chh$t, c(0.025, 0.975)) 
#quantile(boot_overweight_hh_danh$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_dad$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_de$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_goa$t, c(0.025, 0.975)) 
#quantile(boot_overweight_hh_gu$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_ha$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_hp$t, c(0.025, 0.975)) 
#quantile(boot_overweight_hh_jk$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_jh$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_ka$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_ke$t, c(0.025, 0.975)) 
#quantile(boot_overweight_hh_la$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_mp$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_mah$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_man$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_me$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_mi$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_na$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_od$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_pud$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_pun$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_ra$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_si$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_tn$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_te$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_tri$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_up$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_ut$t, c(0.025, 0.975)) 
quantile(boot_overweight_hh_wb$t, c(0.025, 0.975)) 

quantile(boot_obese27_hh_aani$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_ap$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_arp$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_as$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_bi$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_ch$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_chh$t, c(0.025, 0.975)) 
#quantile(boot_obese27_hh_danh$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_dad$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_de$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_goa$t, c(0.025, 0.975)) 
#quantile(boot_obese27_hh_gu$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_ha$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_hp$t, c(0.025, 0.975)) 
#quantile(boot_obese27_hh_jk$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_jh$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_ka$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_ke$t, c(0.025, 0.975)) 
#quantile(boot_obese27_hh_la$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_mp$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_mah$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_man$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_me$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_mi$t, c(0.025, 0.975)) 
 quantile(boot_obese27_hh_na$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_od$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_pud$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_pun$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_ra$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_si$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_tn$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_te$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_tri$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_up$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_ut$t, c(0.025, 0.975)) 
quantile(boot_obese27_hh_wb$t, c(0.025, 0.975)) 

######r ICC by state for clustering on the community level

#DIABETES
diab_psu_aani <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data= AndamanandNicobarIslands)
calc.icc(diab_psu_aani)
diab_psu_ap <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=AndhraPradesh)
calc.icc(diab_psu_ap)
diab_psu_arp <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=ArunachalPradesh)
calc.icc(diab_psu_arp)
diab_psu_as <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Assam)
calc.icc(diab_psu_as)
diab_psu_bi <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Bihar)
calc.icc(diab_psu_bi)
diab_psu_ch <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Chandigarh)
calc.icc(diab_psu_ch)
diab_psu_chh <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Chhattisgarh)
calc.icc(diab_psu_chh)
#diab_psu_danh <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=DadraandNagarHaveli)
#calc.icc(diab_psu_danh)
diab_psu_dad <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=DamanandDiu)
calc.icc(diab_psu_dad)
diab_psu_de <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Delhi)
calc.icc(diab_psu_de)
diab_psu_goa <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Goa)
calc.icc(diab_psu_goa)
#diab_psu_gu <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Gujarat)
#calc.icc(diab_psu_gu)
diab_psu_ha <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Haryana)
calc.icc(diab_psu_ha)
diab_psu_hp <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=HimachalPradesh)
calc.icc(diab_psu_hp)
#diab_psu_jk <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=JammuandKashmir)
#calc.icc(diab_psu_jk)
diab_psu_jh <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Jharkhand)
calc.icc(diab_psu_jh)
diab_psu_ka <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Karnataka)
calc.icc(diab_psu_ka)
diab_psu_ke <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Kerala)
calc.icc(diab_psu_ke)
#diab_psu_la <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Lakshadweep)
#calc.icc(diab_psu_la)
diab_psu_mp <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=MadhyaPradesh)
calc.icc(diab_psu_mp)
diab_psu_mah <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Maharashtra)
calc.icc(diab_psu_mah)
diab_psu_man <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Manipur)
calc.icc(diab_psu_man)
diab_psu_me <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Meghalaya)
calc.icc(diab_psu_me)
diab_psu_mi <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Mizoram)
calc.icc(diab_psu_mi)
diab_psu_na <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Nagaland)
calc.icc(diab_psu_na)
diab_psu_od <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Odisha)
calc.icc(diab_psu_od)
diab_psu_pud <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Puducherry)
calc.icc(diab_psu_pud)
diab_psu_pun <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Punjab)
calc.icc(diab_psu_pun)
diab_psu_ra <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Rajasthan)
calc.icc(diab_psu_ra)
diab_psu_si <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Sikkim)
calc.icc(diab_psu_si)
diab_psu_tn <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=TamilNadu)
calc.icc(diab_psu_tn)
diab_psu_te <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Telangana)
calc.icc(diab_psu_te)
diab_psu_tri <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Tripura)
calc.icc(diab_psu_tri)
diab_psu_up <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=UttarPradesh)
calc.icc(diab_psu_up)
diab_psu_ut <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=Uttarakhand)
calc.icc(diab_psu_ut)
diab_psu_wb <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=WestBengal)
calc.icc(diab_psu_wb)

htn_psu_aani <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data= AndamanandNicobarIslands)
calc.icc(htn_psu_aani)
htn_psu_ap <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=AndhraPradesh)
calc.icc(htn_psu_ap)
htn_psu_arp <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=ArunachalPradesh)
calc.icc(htn_psu_arp)
htn_psu_as <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Assam)
calc.icc(htn_psu_as)
htn_psu_bi <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Bihar)
calc.icc(htn_psu_bi)
htn_psu_ch <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Chandigarh)
calc.icc(htn_psu_ch)
htn_psu_chh <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Chhattisgarh)
calc.icc(htn_psu_chh)
#htn_psu_danh <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=DadraandNagarHaveli)
#calc.icc(htn_psu_danh)
htn_psu_dad <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=DamanandDiu)
calc.icc(htn_psu_dad)
htn_psu_de <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Delhi)
calc.icc(htn_psu_de)
htn_psu_goa <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Goa)
calc.icc(htn_psu_goa)
#htn_psu_gu <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Gujarat)
#calc.icc(htn_psu_gu)
htn_psu_ha <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Haryana)
calc.icc(htn_psu_ha)
htn_psu_hp <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=HimachalPradesh)
calc.icc(htn_psu_hp)
#htn_psu_jk <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=JammuandKashmir)
#calc.icc(htn_psu_jk)
htn_psu_jh <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Jharkhand)
calc.icc(htn_psu_jh)
htn_psu_ka <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Karnataka)
calc.icc(htn_psu_ka)
htn_psu_ke <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Kerala)
calc.icc(htn_psu_ke)
#htn_psu_la <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Lakshadweep)
#calc.icc(htn_psu_la)
htn_psu_mp <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=MadhyaPradesh)
calc.icc(htn_psu_mp)
htn_psu_mah <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Maharashtra)
calc.icc(htn_psu_mah)
htn_psu_man <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Manipur)
calc.icc(htn_psu_man)
htn_psu_me <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Meghalaya)
calc.icc(htn_psu_me)
htn_psu_mi <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Mizoram)
calc.icc(htn_psu_mi)
htn_psu_na <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Nagaland)
calc.icc(htn_psu_na)
htn_psu_od <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Odisha)
calc.icc(htn_psu_od)
htn_psu_pud <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Puducherry)
calc.icc(htn_psu_pud)
htn_psu_pun <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Punjab)
calc.icc(htn_psu_pun)
htn_psu_ra <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Rajasthan)
calc.icc(htn_psu_ra)
htn_psu_si <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Sikkim)
calc.icc(htn_psu_si)
htn_psu_tn <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=TamilNadu)
calc.icc(htn_psu_tn)
htn_psu_te <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Telangana)
calc.icc(htn_psu_te)
htn_psu_tri <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Tripura)
calc.icc(htn_psu_tri)
htn_psu_up <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=UttarPradesh)
calc.icc(htn_psu_up)
htn_psu_ut <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=Uttarakhand)
calc.icc(htn_psu_ut)
htn_psu_wb <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=WestBengal)
calc.icc(htn_psu_wb)


csmoke_psu_aani <- lmer(csmoke ~ 1 + (1|psuid), data= AndamanandNicobarIslands)
calc.icc(csmoke_psu_aani)
csmoke_psu_ap <- lmer(csmoke ~ 1 + (1|psuid), data=AndhraPradesh)
calc.icc(csmoke_psu_ap)
csmoke_psu_arp <- lmer(csmoke ~ 1 + (1|psuid), data=ArunachalPradesh)
calc.icc(csmoke_psu_arp)
csmoke_psu_as <- lmer(csmoke ~ 1 + (1|psuid), data=Assam)
calc.icc(csmoke_psu_as)
csmoke_psu_bi <- lmer(csmoke ~ 1 + (1|psuid), data=Bihar)
calc.icc(csmoke_psu_bi)
csmoke_psu_ch <- lmer(csmoke ~ 1 + (1|psuid), data=Chandigarh)
calc.icc(csmoke_psu_ch)
csmoke_psu_chh <- lmer(csmoke ~ 1 + (1|psuid), data=Chhattisgarh)
calc.icc(csmoke_psu_chh)
#csmoke_psu_danh <- lmer(csmoke ~ 1 + (1|psuid), data=DadraandNagarHaveli)
#calc.icc(csmoke_psu_danh)
csmoke_psu_dad <- lmer(csmoke ~ 1 + (1|psuid), data=DamanandDiu)
calc.icc(csmoke_psu_dad)
csmoke_psu_de <- lmer(csmoke ~ 1 + (1|psuid), data=Delhi)
calc.icc(csmoke_psu_de)
csmoke_psu_goa <- lmer(csmoke ~ 1 + (1|psuid), data=Goa)
calc.icc(csmoke_psu_goa)
#csmoke_psu_gu <- lmer(csmoke ~ 1 + (1|psuid), data=Gujarat)
#calc.icc(csmoke_psu_gu)
csmoke_psu_ha <- lmer(csmoke ~ 1 + (1|psuid), data=Haryana)
calc.icc(csmoke_psu_ha)
csmoke_psu_hp <- lmer(csmoke ~ 1 + (1|psuid), data=HimachalPradesh)
calc.icc(csmoke_psu_hp)
#csmoke_psu_jk <- lmer(csmoke ~ 1 + (1|psuid), data=JammuandKashmir)
#calc.icc(csmoke_psu_jk)
csmoke_psu_jh <- lmer(csmoke ~ 1 + (1|psuid), data=Jharkhand)
calc.icc(csmoke_psu_jh)
csmoke_psu_ka <- lmer(csmoke ~ 1 + (1|psuid), data=Karnataka)
calc.icc(csmoke_psu_ka)
csmoke_psu_ke <- lmer(csmoke ~ 1 + (1|psuid), data=Kerala)
calc.icc(csmoke_psu_ke)
#csmoke_psu_la <- lmer(csmoke ~ 1 + (1|psuid), data=Lakshadweep)
#calc.icc(csmoke_psu_la)
csmoke_psu_mp <- lmer(csmoke ~ 1 + (1|psuid), data=MadhyaPradesh)
calc.icc(csmoke_psu_mp)
csmoke_psu_mah <- lmer(csmoke ~ 1 + (1|psuid), data=Maharashtra)
calc.icc(csmoke_psu_mah)
csmoke_psu_man <- lmer(csmoke ~ 1 + (1|psuid), data=Manipur)
calc.icc(csmoke_psu_man)
csmoke_psu_me <- lmer(csmoke ~ 1 + (1|psuid), data=Meghalaya)
calc.icc(csmoke_psu_me)
csmoke_psu_mi <- lmer(csmoke ~ 1 + (1|psuid), data=Mizoram)
calc.icc(csmoke_psu_mi)
csmoke_psu_na <- lmer(csmoke ~ 1 + (1|psuid), data=Nagaland)
calc.icc(csmoke_psu_na)
csmoke_psu_od <- lmer(csmoke ~ 1 + (1|psuid), data=Odisha)
calc.icc(csmoke_psu_od)
csmoke_psu_pud <- lmer(csmoke ~ 1 + (1|psuid), data=Puducherry)
calc.icc(csmoke_psu_pud)
csmoke_psu_pun <- lmer(csmoke ~ 1 + (1|psuid), data=Punjab)
calc.icc(csmoke_psu_pun)
csmoke_psu_ra <- lmer(csmoke ~ 1 + (1|psuid), data=Rajasthan)
calc.icc(csmoke_psu_ra)
csmoke_psu_si <- lmer(csmoke ~ 1 + (1|psuid), data=Sikkim)
calc.icc(csmoke_psu_si)
csmoke_psu_tn <- lmer(csmoke ~ 1 + (1|psuid), data=TamilNadu)
calc.icc(csmoke_psu_tn)
csmoke_psu_te <- lmer(csmoke ~ 1 + (1|psuid), data=Telangana)
calc.icc(csmoke_psu_te)
csmoke_psu_tri <- lmer(csmoke ~ 1 + (1|psuid), data=Tripura)
calc.icc(csmoke_psu_tri)
csmoke_psu_up <- lmer(csmoke ~ 1 + (1|psuid), data=UttarPradesh)
calc.icc(csmoke_psu_up)
csmoke_psu_ut <- lmer(csmoke ~ 1 + (1|psuid), data=Uttarakhand)
calc.icc(csmoke_psu_ut)
csmoke_psu_wb <- lmer(csmoke ~ 1 + (1|psuid), data=WestBengal)
calc.icc(csmoke_psu_wb)

overweight_psu_aani <- lmer(overweight23 ~ 1 + (1|psuid), data= AndamanandNicobarIslands)
calc.icc(overweight_psu_aani)
overweight_psu_ap <- lmer(overweight23 ~ 1 + (1|psuid), data=AndhraPradesh)
calc.icc(overweight_psu_ap)
overweight_psu_arp <- lmer(overweight23 ~ 1 + (1|psuid), data=ArunachalPradesh)
calc.icc(overweight_psu_arp)
overweight_psu_as <- lmer(overweight23 ~ 1 + (1|psuid), data=Assam)
calc.icc(overweight_psu_as)
overweight_psu_bi <- lmer(overweight23 ~ 1 + (1|psuid), data=Bihar)
calc.icc(overweight_psu_bi)
overweight_psu_ch <- lmer(overweight23 ~ 1 + (1|psuid), data=Chandigarh)
calc.icc(overweight_psu_ch)
overweight_psu_chh <- lmer(overweight23 ~ 1 + (1|psuid), data=Chhattisgarh)
calc.icc(overweight_psu_chh)
#overweight_psu_danh <- lmer(overweight23 ~ 1 + (1|psuid), data=DadraandNagarHaveli)
#calc.icc(overweight_psu_danh)
overweight_psu_dad <- lmer(overweight23 ~ 1 + (1|psuid), data=DamanandDiu)
calc.icc(overweight_psu_dad)
overweight_psu_de <- lmer(overweight23 ~ 1 + (1|psuid), data=Delhi)
calc.icc(overweight_psu_de)
overweight_psu_goa <- lmer(overweight23 ~ 1 + (1|psuid), data=Goa)
calc.icc(overweight_psu_goa)
#overweight_psu_gu <- lmer(overweight23 ~ 1 + (1|psuid), data=Gujarat)
#calc.icc(overweight_psu_gu)
overweight_psu_ha <- lmer(overweight23 ~ 1 + (1|psuid), data=Haryana)
calc.icc(overweight_psu_ha)
overweight_psu_hp <- lmer(overweight23 ~ 1 + (1|psuid), data=HimachalPradesh)
calc.icc(overweight_psu_hp)
#overweight_psu_jk <- lmer(overweight23 ~ 1 + (1|psuid), data=JammuandKashmir)
#calc.icc(overweight_psu_jk)
overweight_psu_jh <- lmer(overweight23 ~ 1 + (1|psuid), data=Jharkhand)
calc.icc(overweight_psu_jh)
overweight_psu_ka <- lmer(overweight23 ~ 1 + (1|psuid), data=Karnataka)
calc.icc(overweight_psu_ka)
overweight_psu_ke <- lmer(overweight23 ~ 1 + (1|psuid), data=Kerala)
calc.icc(overweight_psu_ke)
#overweight_psu_la <- lmer(overweight23 ~ 1 + (1|psuid), data=Lakshadweep)
#calc.icc(overweight_psu_la)
overweight_psu_mp <- lmer(overweight23 ~ 1 + (1|psuid), data=MadhyaPradesh)
calc.icc(overweight_psu_mp)
overweight_psu_mah <- lmer(overweight23 ~ 1 + (1|psuid), data=Maharashtra)
calc.icc(overweight_psu_mah)
overweight_psu_man <- lmer(overweight23 ~ 1 + (1|psuid), data=Manipur)
calc.icc(overweight_psu_man)
overweight_psu_me <- lmer(overweight23 ~ 1 + (1|psuid), data=Meghalaya)
calc.icc(overweight_psu_me)
overweight_psu_mi <- lmer(overweight23 ~ 1 + (1|psuid), data=Mizoram)
calc.icc(overweight_psu_mi)
overweight_psu_na <- lmer(overweight23 ~ 1 + (1|psuid), data=Nagaland)
calc.icc(overweight_psu_na)
overweight_psu_od <- lmer(overweight23 ~ 1 + (1|psuid), data=Odisha)
calc.icc(overweight_psu_od)
overweight_psu_pud <- lmer(overweight23 ~ 1 + (1|psuid), data=Puducherry)
calc.icc(overweight_psu_pud)
overweight_psu_pun <- lmer(overweight23 ~ 1 + (1|psuid), data=Punjab)
calc.icc(overweight_psu_pun)
overweight_psu_ra <- lmer(overweight23 ~ 1 + (1|psuid), data=Rajasthan)
calc.icc(overweight_psu_ra)
overweight_psu_si <- lmer(overweight23 ~ 1 + (1|psuid), data=Sikkim)
calc.icc(overweight_psu_si)
overweight_psu_tn <- lmer(overweight23 ~ 1 + (1|psuid), data=TamilNadu)
calc.icc(overweight_psu_tn)
overweight_psu_te <- lmer(overweight23 ~ 1 + (1|psuid), data=Telangana)
calc.icc(overweight_psu_te)
overweight_psu_tri <- lmer(overweight23 ~ 1 + (1|psuid), data=Tripura)
calc.icc(overweight_psu_tri)
overweight_psu_up <- lmer(overweight23 ~ 1 + (1|psuid), data=UttarPradesh)
calc.icc(overweight_psu_up)
overweight_psu_ut <- lmer(overweight23 ~ 1 + (1|psuid), data=Uttarakhand)
calc.icc(overweight_psu_ut)
overweight_psu_wb <- lmer(overweight23 ~ 1 + (1|psuid), data=WestBengal)
calc.icc(overweight_psu_wb)


obese27_psu_aani <- lmer(obese27 ~ 1 + (1|psuid), data= AndamanandNicobarIslands)
calc.icc(obese27_psu_aani)
obese27_psu_ap <- lmer(obese27 ~ 1 + (1|psuid), data=AndhraPradesh)
calc.icc(obese27_psu_ap)
obese27_psu_arp <- lmer(obese27 ~ 1 + (1|psuid), data=ArunachalPradesh)
calc.icc(obese27_psu_arp)
obese27_psu_as <- lmer(obese27 ~ 1 + (1|psuid), data=Assam)
calc.icc(obese27_psu_as)
obese27_psu_bi <- lmer(obese27 ~ 1 + (1|psuid), data=Bihar)
calc.icc(obese27_psu_bi)
obese27_psu_ch <- lmer(obese27 ~ 1 + (1|psuid), data=Chandigarh)
calc.icc(obese27_psu_ch)
obese27_psu_chh <- lmer(obese27 ~ 1 + (1|psuid), data=Chhattisgarh)
calc.icc(obese27_psu_chh)
#obese27_psu_danh <- lmer(obese27 ~ 1 + (1|psuid), data=DadraandNagarHaveli)
#calc.icc(obese27_psu_danh)
obese27_psu_dad <- lmer(obese27 ~ 1 + (1|psuid), data=DamanandDiu)
calc.icc(obese27_psu_dad)
obese27_psu_de <- lmer(obese27 ~ 1 + (1|psuid), data=Delhi)
calc.icc(obese27_psu_de)
obese27_psu_goa <- lmer(obese27 ~ 1 + (1|psuid), data=Goa)
calc.icc(obese27_psu_goa)
#obese27_psu_gu <- lmer(obese27 ~ 1 + (1|psuid), data=Gujarat)
#calc.icc(obese27_psu_gu)
obese27_psu_ha <- lmer(obese27 ~ 1 + (1|psuid), data=Haryana)
calc.icc(obese27_psu_ha)
obese27_psu_hp <- lmer(obese27 ~ 1 + (1|psuid), data=HimachalPradesh)
calc.icc(obese27_psu_hp)
#obese27_psu_jk <- lmer(obese27 ~ 1 + (1|psuid), data=JammuandKashmir)
#calc.icc(obese27_psu_jk)
obese27_psu_jh <- lmer(obese27 ~ 1 + (1|psuid), data=Jharkhand)
calc.icc(obese27_psu_jh)
obese27_psu_ka <- lmer(obese27 ~ 1 + (1|psuid), data=Karnataka)
calc.icc(obese27_psu_ka)
obese27_psu_ke <- lmer(obese27 ~ 1 + (1|psuid), data=Kerala)
calc.icc(obese27_psu_ke)
#obese27_psu_la <- lmer(obese27 ~ 1 + (1|psuid), data=Lakshadweep)
#calc.icc(obese27_psu_la)
obese27_psu_mp <- lmer(obese27 ~ 1 + (1|psuid), data=MadhyaPradesh)
calc.icc(obese27_psu_mp)
obese27_psu_mah <- lmer(obese27 ~ 1 + (1|psuid), data=Maharashtra)
calc.icc(obese27_psu_mah)
obese27_psu_man <- lmer(obese27 ~ 1 + (1|psuid), data=Manipur)
calc.icc(obese27_psu_man)
obese27_psu_me <- lmer(obese27 ~ 1 + (1|psuid), data=Meghalaya)
calc.icc(obese27_psu_me)
obese27_psu_mi <- lmer(obese27 ~ 1 + (1|psuid), data=Mizoram)
calc.icc(obese27_psu_mi)
obese27_psu_na <- lmer(obese27 ~ 1 + (1|psuid), data=Nagaland)
calc.icc(obese27_psu_na)
obese27_psu_od <- lmer(obese27 ~ 1 + (1|psuid), data=Odisha)
calc.icc(obese27_psu_od)
obese27_psu_pud <- lmer(obese27 ~ 1 + (1|psuid), data=Puducherry)
calc.icc(obese27_psu_pud)
obese27_psu_pun <- lmer(obese27 ~ 1 + (1|psuid), data=Punjab)
calc.icc(obese27_psu_pun)
obese27_psu_ra <- lmer(obese27 ~ 1 + (1|psuid), data=Rajasthan)
calc.icc(obese27_psu_ra)
obese27_psu_si <- lmer(obese27 ~ 1 + (1|psuid), data=Sikkim)
calc.icc(obese27_psu_si)
obese27_psu_tn <- lmer(obese27 ~ 1 + (1|psuid), data=TamilNadu)
calc.icc(obese27_psu_tn)
obese27_psu_te <- lmer(obese27 ~ 1 + (1|psuid), data=Telangana)
calc.icc(obese27_psu_te)
obese27_psu_tri <- lmer(obese27 ~ 1 + (1|psuid), data=Tripura)
calc.icc(obese27_psu_tri)
obese27_psu_up <- lmer(obese27 ~ 1 + (1|psuid), data=UttarPradesh)
calc.icc(obese27_psu_up)
obese27_psu_ut <- lmer(obese27 ~ 1 + (1|psuid), data=Uttarakhand)
calc.icc(obese27_psu_ut)
obese27_psu_wb <- lmer(obese27 ~ 1 + (1|psuid), data=WestBengal)
calc.icc(obese27_psu_wb)



############calculate bootstraps for state iccs for community clustering########

boot_diab_psu_aani <- bootMer(diab_psu_aani, calc.icc, nsim=500)
boot_diab_psu_ap <- bootMer(diab_psu_ap, calc.icc, nsim=500)
boot_diab_psu_arp <- bootMer(diab_psu_arp, calc.icc, nsim=500)
boot_diab_psu_as <- bootMer(diab_psu_as, calc.icc, nsim=500)
boot_diab_psu_bi <- bootMer(diab_psu_bi, calc.icc, nsim=500)
boot_diab_psu_ch <- bootMer(diab_psu_ch, calc.icc, nsim=500)
boot_diab_psu_chh <- bootMer(diab_psu_chh, calc.icc, nsim=500)
#boot_diab_psu_danh <- bootMer(diab_psu_danh, calc.icc, nsim=500)
boot_diab_psu_dad <- bootMer(diab_psu_dad, calc.icc, nsim=500)
boot_diab_psu_de <- bootMer(diab_psu_de, calc.icc, nsim=500) 
boot_diab_psu_goa <- bootMer(diab_psu_goa, calc.icc, nsim=500)
#boot_diab_psu_gu <- bootMer(diab_psu_gu, calc.icc, nsim=500)
boot_diab_psu_ha <- bootMer(diab_psu_ha, calc.icc, nsim=500)
boot_diab_psu_hp <- bootMer(diab_psu_hp, calc.icc, nsim=500)
#boot_diab_psu_jk <- bootMer(diab_psu_jk, calc.icc, nsim=500)
boot_diab_psu_jh <- bootMer(diab_psu_jh, calc.icc, nsim=500)
boot_diab_psu_ka <- bootMer(diab_psu_ka, calc.icc, nsim=500)
boot_diab_psu_ke <- bootMer(diab_psu_ke, calc.icc, nsim=500)
#boot_diab_psu_la <- bootMer(diab_psu_la, calc.icc, nsim=500)
boot_diab_psu_mp <- bootMer(diab_psu_mp, calc.icc, nsim=500)
boot_diab_psu_mah <- bootMer(diab_psu_mah, calc.icc, nsim=500)
boot_diab_psu_man <- bootMer(diab_psu_man, calc.icc, nsim=500)
boot_diab_psu_me <- bootMer(diab_psu_me, calc.icc, nsim=500)
boot_diab_psu_mi <- bootMer(diab_psu_mi, calc.icc, nsim=500) 
boot_diab_psu_na <- bootMer(diab_psu_na, calc.icc, nsim=500)
boot_diab_psu_od <- bootMer(diab_psu_od, calc.icc, nsim=500)
boot_diab_psu_pud <- bootMer(diab_psu_pud, calc.icc, nsim=500)
boot_diab_psu_pun <- bootMer(diab_psu_pun, calc.icc, nsim=500)
boot_diab_psu_ra <- bootMer(diab_psu_ra, calc.icc, nsim=500)
boot_diab_psu_si <- bootMer(diab_psu_si, calc.icc, nsim=500)
boot_diab_psu_tn <- bootMer(diab_psu_tn, calc.icc, nsim=500)
boot_diab_psu_te <- bootMer(diab_psu_te, calc.icc, nsim=500)
boot_diab_psu_tri <- bootMer(diab_psu_tri, calc.icc, nsim=500)
boot_diab_psu_up <- bootMer(diab_psu_up, calc.icc, nsim=500)
boot_diab_psu_ut <- bootMer(diab_psu_ut, calc.icc, nsim=500)
boot_diab_psu_wb <- bootMer(diab_psu_wb, calc.icc, nsim=500)

boot_htn_psu_aani <- bootMer(htn_psu_aani, calc.icc, nsim=500)
boot_htn_psu_ap <- bootMer(htn_psu_ap, calc.icc, nsim=500)
boot_htn_psu_arp <- bootMer(htn_psu_arp, calc.icc, nsim=500)
boot_htn_psu_as <- bootMer(htn_psu_as, calc.icc, nsim=500)
boot_htn_psu_bi <- bootMer(htn_psu_bi, calc.icc, nsim=500)
boot_htn_psu_ch <- bootMer(htn_psu_ch, calc.icc, nsim=500)
boot_htn_psu_chh <- bootMer(htn_psu_chh, calc.icc, nsim=500)
#boot_htn_psu_danh <- bootMer(htn_psu_danh, calc.icc, nsim=500)
boot_htn_psu_dad <- bootMer(htn_psu_dad, calc.icc, nsim=500)
boot_htn_psu_de <- bootMer(htn_psu_de, calc.icc, nsim=500)
boot_htn_psu_goa <- bootMer(htn_psu_goa, calc.icc, nsim=500)
#boot_htn_psu_gu <- bootMer(htn_psu_gu, calc.icc, nsim=500)
boot_htn_psu_ha <- bootMer(htn_psu_ha, calc.icc, nsim=500)
boot_htn_psu_hp <- bootMer(htn_psu_hp, calc.icc, nsim=500)
#boot_htn_psu_jk <- bootMer(htn_psu_jk, calc.icc, nsim=500)
boot_htn_psu_jh <- bootMer(htn_psu_jh, calc.icc, nsim=500)
boot_htn_psu_ka <- bootMer(htn_psu_ka, calc.icc, nsim=500)
boot_htn_psu_ke <- bootMer(htn_psu_ke, calc.icc, nsim=500)
#boot_htn_psu_la <- bootMer(htn_psu_la, calc.icc, nsim=500)
boot_htn_psu_mp <- bootMer(htn_psu_mp, calc.icc, nsim=500)
boot_htn_psu_mah <- bootMer(htn_psu_mah, calc.icc, nsim=500)
boot_htn_psu_man <- bootMer(htn_psu_man, calc.icc, nsim=500)
boot_htn_psu_me <- bootMer(htn_psu_me, calc.icc, nsim=500)
boot_htn_psu_mi <- bootMer(htn_psu_mi, calc.icc, nsim=500)
boot_htn_psu_na <- bootMer(htn_psu_na, calc.icc, nsim=500)
boot_htn_psu_od <- bootMer(htn_psu_od, calc.icc, nsim=500)
boot_htn_psu_pud <- bootMer(htn_psu_pud, calc.icc, nsim=500)
boot_htn_psu_pun <- bootMer(htn_psu_pun, calc.icc, nsim=500)
boot_htn_psu_ra <- bootMer(htn_psu_ra, calc.icc, nsim=500)
boot_htn_psu_si <- bootMer(htn_psu_si, calc.icc, nsim=500)
boot_htn_psu_tn <- bootMer(htn_psu_tn, calc.icc, nsim=500)
boot_htn_psu_te <- bootMer(htn_psu_te, calc.icc, nsim=500)
boot_htn_psu_tri <- bootMer(htn_psu_tri, calc.icc, nsim=500)
boot_htn_psu_up <- bootMer(htn_psu_up, calc.icc, nsim=500)
boot_htn_psu_ut <- bootMer(htn_psu_ut, calc.icc, nsim=500)
boot_htn_psu_wb <- bootMer(htn_psu_wb, calc.icc, nsim=500)

boot_csmoke_psu_aani <- bootMer(csmoke_psu_aani, calc.icc, nsim=500)
boot_csmoke_psu_ap <- bootMer(csmoke_psu_ap, calc.icc, nsim=500)
boot_csmoke_psu_arp <- bootMer(csmoke_psu_arp, calc.icc, nsim=500)
boot_csmoke_psu_as <- bootMer(csmoke_psu_as, calc.icc, nsim=500)
boot_csmoke_psu_bi <- bootMer(csmoke_psu_bi, calc.icc, nsim=500)
boot_csmoke_psu_ch <- bootMer(csmoke_psu_ch, calc.icc, nsim=500)
boot_csmoke_psu_chh <- bootMer(csmoke_psu_chh, calc.icc, nsim=500)
#boot_csmoke_psu_danh <- bootMer(csmoke_psu_danh, calc.icc, nsim=500)
boot_csmoke_psu_dad <- bootMer(csmoke_psu_dad, calc.icc, nsim=500)
boot_csmoke_psu_de <- bootMer(csmoke_psu_de, calc.icc, nsim=500)
boot_csmoke_psu_goa <- bootMer(csmoke_psu_goa, calc.icc, nsim=500)
#boot_csmoke_psu_gu <- bootMer(csmoke_psu_gu, calc.icc, nsim=500)
boot_csmoke_psu_ha <- bootMer(csmoke_psu_ha, calc.icc, nsim=500)
boot_csmoke_psu_hp <- bootMer(csmoke_psu_hp, calc.icc, nsim=500)
#boot_csmoke_psu_jk <- bootMer(csmoke_psu_jk, calc.icc, nsim=500)
boot_csmoke_psu_jh <- bootMer(csmoke_psu_jh, calc.icc, nsim=500)
boot_csmoke_psu_ka <- bootMer(csmoke_psu_ka, calc.icc, nsim=500)
boot_csmoke_psu_ke <- bootMer(csmoke_psu_ke, calc.icc, nsim=500)
#boot_csmoke_psu_la <- bootMer(csmoke_psu_la, calc.icc, nsim=500)
boot_csmoke_psu_mp <- bootMer(csmoke_psu_mp, calc.icc, nsim=500)
boot_csmoke_psu_mah <- bootMer(csmoke_psu_mah, calc.icc, nsim=500)
boot_csmoke_psu_man <- bootMer(csmoke_psu_man, calc.icc, nsim=500)
boot_csmoke_psu_me <- bootMer(csmoke_psu_me, calc.icc, nsim=500)
boot_csmoke_psu_mi <- bootMer(csmoke_psu_mi, calc.icc, nsim=500)
boot_csmoke_psu_na <- bootMer(csmoke_psu_na, calc.icc, nsim=500)
boot_csmoke_psu_od <- bootMer(csmoke_psu_od, calc.icc, nsim=500)
boot_csmoke_psu_pud <- bootMer(csmoke_psu_pud, calc.icc, nsim=500)
boot_csmoke_psu_pun <- bootMer(csmoke_psu_pun, calc.icc, nsim=500)
boot_csmoke_psu_ra <- bootMer(csmoke_psu_ra, calc.icc, nsim=500)
boot_csmoke_psu_si <- bootMer(csmoke_psu_si, calc.icc, nsim=500)
boot_csmoke_psu_tn <- bootMer(csmoke_psu_tn, calc.icc, nsim=500)
boot_csmoke_psu_te <- bootMer(csmoke_psu_te, calc.icc, nsim=500)
boot_csmoke_psu_tri <- bootMer(csmoke_psu_tri, calc.icc, nsim=500)
boot_csmoke_psu_up <- bootMer(csmoke_psu_up, calc.icc, nsim=500)
boot_csmoke_psu_ut <- bootMer(csmoke_psu_ut, calc.icc, nsim=500)
boot_csmoke_psu_wb <- bootMer(csmoke_psu_wb, calc.icc, nsim=500)

boot_overweight_psu_aani <- bootMer(overweight_psu_aani, calc.icc, nsim=500)
boot_overweight_psu_ap <- bootMer(overweight_psu_ap, calc.icc, nsim=500)
boot_overweight_psu_arp <- bootMer(overweight_psu_arp, calc.icc, nsim=500)
boot_overweight_psu_as <- bootMer(overweight_psu_as, calc.icc, nsim=500)
boot_overweight_psu_bi <- bootMer(overweight_psu_bi, calc.icc, nsim=500)
boot_overweight_psu_ch <- bootMer(overweight_psu_ch, calc.icc, nsim=500)
boot_overweight_psu_chh <- bootMer(overweight_psu_chh, calc.icc, nsim=500)
#boot_overweight_psu_danh <- bootMer(overweight_psu_danh, calc.icc, nsim=500)
boot_overweight_psu_dad <- bootMer(overweight_psu_dad, calc.icc, nsim=500)
boot_overweight_psu_de <- bootMer(overweight_psu_de, calc.icc, nsim=500)
boot_overweight_psu_goa <- bootMer(overweight_psu_goa, calc.icc, nsim=500)
#boot_overweight_psu_gu <- bootMer(overweight_psu_gu, calc.icc, nsim=500)
boot_overweight_psu_ha <- bootMer(overweight_psu_ha, calc.icc, nsim=500)
boot_overweight_psu_hp <- bootMer(overweight_psu_hp, calc.icc, nsim=500)
#boot_overweight_psu_jk <- bootMer(overweight_psu_jk, calc.icc, nsim=500)
boot_overweight_psu_jh <- bootMer(overweight_psu_jh, calc.icc, nsim=500)
boot_overweight_psu_ka <- bootMer(overweight_psu_ka, calc.icc, nsim=500)
boot_overweight_psu_ke <- bootMer(overweight_psu_ke, calc.icc, nsim=500)
#boot_overweight_psu_la <- bootMer(overweight_psu_la, calc.icc, nsim=500)
boot_overweight_psu_mp <- bootMer(overweight_psu_mp, calc.icc, nsim=500)
boot_overweight_psu_mah <- bootMer(overweight_psu_mah, calc.icc, nsim=500)
boot_overweight_psu_man <- bootMer(overweight_psu_man, calc.icc, nsim=500)
boot_overweight_psu_me <- bootMer(overweight_psu_me, calc.icc, nsim=500)
boot_overweight_psu_mi <- bootMer(overweight_psu_mi, calc.icc, nsim=500)
boot_overweight_psu_na <- bootMer(overweight_psu_na, calc.icc, nsim=500)
boot_overweight_psu_od <- bootMer(overweight_psu_od, calc.icc, nsim=500)
boot_overweight_psu_pud <- bootMer(overweight_psu_pud, calc.icc, nsim=500)
boot_overweight_psu_pun <- bootMer(overweight_psu_pun, calc.icc, nsim=500)
boot_overweight_psu_ra <- bootMer(overweight_psu_ra, calc.icc, nsim=500)
boot_overweight_psu_si <- bootMer(overweight_psu_si, calc.icc, nsim=500)
boot_overweight_psu_tn <- bootMer(overweight_psu_tn, calc.icc, nsim=500)
boot_overweight_psu_te <- bootMer(overweight_psu_te, calc.icc, nsim=500)
boot_overweight_psu_tri <- bootMer(overweight_psu_tri, calc.icc, nsim=500)
boot_overweight_psu_up <- bootMer(overweight_psu_up, calc.icc, nsim=500)
boot_overweight_psu_ut <- bootMer(overweight_psu_ut, calc.icc, nsim=500)
boot_overweight_psu_wb <- bootMer(overweight_psu_wb, calc.icc, nsim=500)

boot_obese27_psu_aani <- bootMer(obese27_psu_aani, calc.icc, nsim=500)
boot_obese27_psu_ap <- bootMer(obese27_psu_ap, calc.icc, nsim=500)
boot_obese27_psu_arp <- bootMer(obese27_psu_arp, calc.icc, nsim=500)
boot_obese27_psu_as <- bootMer(obese27_psu_as, calc.icc, nsim=500)
boot_obese27_psu_bi <- bootMer(obese27_psu_bi, calc.icc, nsim=500)
boot_obese27_psu_ch <- bootMer(obese27_psu_ch, calc.icc, nsim=500)
boot_obese27_psu_chh <- bootMer(obese27_psu_chh, calc.icc, nsim=500)
#boot_obese27_psu_danh <- bootMer(obese27_psu_danh, calc.icc, nsim=500)
boot_obese27_psu_dad <- bootMer(obese27_psu_dad, calc.icc, nsim=500)
boot_obese27_psu_de <- bootMer(obese27_psu_de, calc.icc, nsim=500)
boot_obese27_psu_goa <- bootMer(obese27_psu_goa, calc.icc, nsim=500)
#boot_obese27_psu_gu <- bootMer(obese27_psu_gu, calc.icc, nsim=500)
boot_obese27_psu_ha <- bootMer(obese27_psu_ha, calc.icc, nsim=500)
boot_obese27_psu_hp <- bootMer(obese27_psu_hp, calc.icc, nsim=500)
#boot_obese27_psu_jk <- bootMer(obese27_psu_jk ,calc.icc, nsim=500)
boot_obese27_psu_jh <- bootMer(obese27_psu_jh, calc.icc, nsim=500)
boot_obese27_psu_ka <- bootMer(obese27_psu_ka, calc.icc, nsim=500)
boot_obese27_psu_ke <- bootMer(obese27_psu_ke, calc.icc, nsim=500)
#boot_obese27_psu_la <- bootMer(obese27_psu_la, calc.icc, nsim=500)
boot_obese27_psu_mp <- bootMer(obese27_psu_mp, calc.icc, nsim=500)
boot_obese27_psu_mah <- bootMer(obese27_psu_mah, calc.icc, nsim=500)
boot_obese27_psu_man <- bootMer(obese27_psu_man, calc.icc, nsim=500)
boot_obese27_psu_me <- bootMer(obese27_psu_me, calc.icc, nsim=500)
boot_obese27_psu_mi <- bootMer(obese27_psu_mi, calc.icc, nsim=500)
boot_obese27_psu_na <- bootMer(obese27_psu_na, calc.icc, nsim=500)
boot_obese27_psu_od <- bootMer(obese27_psu_od, calc.icc, nsim=500)
boot_obese27_psu_pud <- bootMer(obese27_psu_pud, calc.icc, nsim=500)
boot_obese27_psu_pun <- bootMer(obese27_psu_pun, calc.icc, nsim=500)
boot_obese27_psu_ra <- bootMer(obese27_psu_ra, calc.icc, nsim=500)
boot_obese27_psu_si <- bootMer(obese27_psu_si, calc.icc, nsim=500)
boot_obese27_psu_tn <- bootMer(obese27_psu_tn, calc.icc, nsim=500)
boot_obese27_psu_te <- bootMer(obese27_psu_te, calc.icc, nsim=500)
boot_obese27_psu_tri <- bootMer(obese27_psu_tri, calc.icc, nsim=500)
boot_obese27_psu_up <- bootMer(obese27_psu_up, calc.icc, nsim=500)
boot_obese27_psu_ut <- bootMer(obese27_psu_ut, calc.icc, nsim=500)
boot_obese27_psu_wb <- bootMer(obese27_psu_wb, calc.icc, nsim=500)

############################################

quantile(boot_diab_psu_aani$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_ap$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_arp$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_as$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_bi$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_ch$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_chh$t, c(0.025, 0.975)) 
#quantile(boot_diab_psu_danh$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_dad$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_de$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_goa$t, c(0.025, 0.975)) 
#quantile(boot_diab_psu_gu$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_ha$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_hp$t, c(0.025, 0.975)) 
#quantile(boot_diab_psu_jk$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_jh$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_ka$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_ke$t, c(0.025, 0.975)) 
#quantile(boot_diab_psu_la$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_mp$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_mah$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_man$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_me$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_mi$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_na$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_od$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_pud$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_pun$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_ra$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_si$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_tn$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_te$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_tri$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_up$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_ut$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_wb$t, c(0.025, 0.975)) 

quantile(boot_htn_psu_aani$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_ap$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_arp$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_as$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_bi$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_ch$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_chh$t, c(0.025, 0.975)) 
#quantile(boot_htn_psu_danh$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_dad$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_de$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_goa$t, c(0.025, 0.975)) 
#quantile(boot_htn_psu_gu$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_ha$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_hp$t, c(0.025, 0.975)) 
#quantile(boot_htn_psu_jk$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_jh$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_ka$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_ke$t, c(0.025, 0.975)) 
#quantile(boot_htn_psu_la$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_mp$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_mah$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_man$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_me$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_mi$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_na$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_od$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_pud$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_pun$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_ra$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_si$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_tn$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_te$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_tri$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_up$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_ut$t, c(0.025, 0.975)) 
quantile(boot_htn_psu_wb$t, c(0.025, 0.975)) 

quantile(boot_csmoke_psu_aani$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_ap$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_arp$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_as$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_bi$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_ch$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_chh$t, c(0.025, 0.975)) 
#quantile(boot_csmoke_psu_danh$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_dad$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_de$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_goa$t, c(0.025, 0.975)) 
#quantile(boot_csmoke_psu_gu$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_ha$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_hp$t, c(0.025, 0.975)) 
#quantile(boot_csmoke_psu_jk$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_jh$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_ka$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_ke$t, c(0.025, 0.975)) 
#quantile(boot_csmoke_psu_la$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_mp$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_mah$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_man$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_me$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_mi$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_na$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_od$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_pud$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_pun$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_ra$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_si$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_tn$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_te$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_tri$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_up$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_ut$t, c(0.025, 0.975)) 
quantile(boot_csmoke_psu_wb$t, c(0.025, 0.975)) 

quantile(boot_overweight_psu_aani$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_ap$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_arp$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_as$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_bi$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_ch$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_chh$t, c(0.025, 0.975)) 
#quantile(boot_overweight_psu_danh$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_dad$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_de$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_goa$t, c(0.025, 0.975)) 
#quantile(boot_overweight_psu_gu$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_ha$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_hp$t, c(0.025, 0.975)) 
#quantile(boot_overweight_psu_jk$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_jh$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_ka$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_ke$t, c(0.025, 0.975)) 
#quantile(boot_overweight_psu_la$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_mp$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_mah$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_man$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_me$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_mi$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_na$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_od$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_pud$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_pun$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_ra$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_si$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_tn$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_te$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_tri$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_up$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_ut$t, c(0.025, 0.975)) 
quantile(boot_overweight_psu_wb$t, c(0.025, 0.975)) 

quantile(boot_obese27_psu_aani$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_ap$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_arp$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_as$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_bi$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_ch$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_chh$t, c(0.025, 0.975)) 
#quantile(boot_obese27_psu_danh$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_dad$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_de$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_goa$t, c(0.025, 0.975)) 
#quantile(boot_obese27_psu_gu$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_ha$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_hp$t, c(0.025, 0.975)) 
#quantile(boot_obese27_psu_jk$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_jh$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_ka$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_ke$t, c(0.025, 0.975)) 
#quantile(boot_obese27_psu_la$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_mp$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_mah$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_man$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_me$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_mi$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_na$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_od$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_pud$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_pun$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_ra$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_si$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_tn$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_te$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_tri$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_up$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_ut$t, c(0.025, 0.975)) 
quantile(boot_obese27_psu_wb$t, c(0.025, 0.975))
```


```{Figure 2: household ICC by district~ median district wealth index}

# calculates lmer model for each d_id subset for ICC on hh_id level

india$ex_diab_narrow_ind<-as.numeric(levels(india$ex_diab_narrow_ind))[india$ex_diab_narrow_ind]
india$ex_htn_narrow_ind<-as.numeric(levels(india$ex_htn_narrow_ind))[india$ex_htn_narrow_ind]
india$csmoke<-as.numeric(levels(india$csmoke))[india$csmoke]

#calculate median of household wealth index
library(spatstat)

india <- india %>% 
  mutate(urban_dbl = as.numeric(urban)-1)

india <- india %>% 
  mutate(urban_lab = as.factor(ifelse(is.na(rural) == TRUE, NA, 
                                  ifelse(rural == 1, "Rural", "Urban"))))

distmean.dat <- india %>%
  group_by(d_id) %>% 
  mutate(prop_urban = 100*weighted.mean(urban_dbl, sworld_weight_india, na.rm=TRUE))%>% 
  ungroup() %>% 
  group_by(d_id, urban_lab) %>%
  mutate(medianai = spatstat::weighted.median(ai_NAT_rurb, sworld_weight_india, na.rm=TRUE)) %>%
  filter(row_number()==1) %>%
  ungroup() %>% 
  #mutate(q_medianai = ntile(medianai, 5)) %>% 
  dplyr::select( medianai, zone,d_id, urban_lab, prop_urban)
#use ai_NAT_rurb because the figures are created separately for rural/urban

#before calculating ICCs by district, district with <50 obs filtered out (ICCs would be 0/1/NA)
india$d_id<-as.numeric(levels(india$d_id))[india$d_id]
india_50d <- india %>%
  group_by(d_id) %>%
  mutate(n_district=n()) %>%
  filter (n_district >=50) #1 district filtered (d_id=1823)

#then reconvert to vector for loop
india_50d$d_id <- as.factor(india_50d$d_id)

#add columns with amount of diabetics/hypertensive people etc per district

library(dplyr)

india_50d_diab <-india_50d %>% 
  group_by(d_id) %>% 
  summarise(n_diab = sum(ex_diab_narrow_ind))

india_50d_htn <-india_50d %>% 
  group_by(d_id) %>% 
  summarise(n_htn = sum(ex_htn_narrow_ind))
   
india_50d_smok <-india_50d %>% 
  group_by(d_id) %>% 
  summarise(n_smok = sum(csmoke))

india_50d_ovw <-india_50d %>% 
  group_by(d_id) %>% 
  summarise(n_ovw = sum(overweight23))

india_50d_obese <-india_50d %>% 
  group_by(d_id) %>% 
  summarise(n_obese = sum(obese27))

india_50d<- left_join (india_50d, india_50d_diab)
india_50d<- left_join (india_50d, india_50d_htn)
india_50d<- left_join (india_50d, india_50d_smok)
india_50d<- left_join (india_50d, india_50d_ovw)
india_50d<- left_join (india_50d, india_50d_obese)

#separate dataset into rural and urban

rural <- filter(india_50d, urban==0)
urban <- filter(india_50d, urban==1)

#rural diabetes

# filter out districts with less than 20 diabetics (ICCs showed to be 0/1)

library(dplyr)

rural$d_id<-as.numeric(levels(rural$d_id))[rural$d_id]

rural_diab <- rural %>%
  filter(n_diab >=20)

rural_diab$d_id <- as.factor(rural_diab$d_id)

for (i in levels(rural_diab$d_id)){
  print(paste( i))
  df1 = subset(rural_diab, d_id == i)
  model <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data= df1 )
  icc<-calc.icc(model)
  print(icc)
} 
#gives a row with district followed by a row with the ICC, formatted in excel


#rural htn
rural$d_id<-as.numeric(levels(rural$d_id))[rural$d_id]

rural_htn <- rural %>%
  filter(n_htn >=20)

rural_htn$d_id <- as.factor(rural_htn$d_id)

for (i in levels(rural_htn$d_id)){
  print(paste( i))
  df1 = subset(rural_htn, d_id == i)
  model <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data= df1 )
  icc<-calc.icc(model)
  print(icc)
} 

#rural smoking

rural_smok <- rural %>%
  filter(n_htn >=20)

rural_smok$d_id <- as.factor(rural_smok$d_id)    
  
for (i in levels(rural_smok$d_id)){
  print(paste( i))
  df1 = subset(rural_smok, d_id == i)
  model <- lmer(csmoke ~ 1 + (1|hh_id), data= df1 )
  icc<-calc.icc(model)
  print(icc)
} 

#rural overweight

rural_ovw <- rural %>%
  filter(n_ovw >=20)

rural_ovw$d_id <- as.factor(rural_ovw$d_id) 

for (i in levels(rural_ovw$d_id)){
  print(paste( i))
  df1 = subset(rural:ovw, d_id == i)
  model <- lmer(overweight23~ 1 + (1|hh_id), data= df1 )
  icc<-calc.icc(model)
  print(icc)
} 

#rural obese

rural_obese <- rural %>%
  filter(n_obese >=20)

rural_obese$d_id <- as.factor(rural_obese$d_id) 

for (i in levels(rural_obese$d_id)){
  print(paste( i))
  df1 = subset(rural_obese, d_id == i)
  model <- lmer(obese27 ~ 1 + (1|hh_id), data= df1 )
  icc<-calc.icc(model)
  print(icc)
} 

######urban########

#urban/diabetes

urban$d_id<-as.numeric(levels(urban$d_id))[urban$d_id]

urban_diab <- urban %>%
  filter(n_diab >=20)

urban_diab$d_id <- as.factor(urban_diab$d_id)

for (i in levels(urban_diab$d_id)){
  print(paste( i))
  df1 = subset(urban_diab, d_id == i)
  model <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data= df1 )
  icc<-calc.icc(model)
  print(icc)
} 

#urban/hypertension

urban_htn <- urban %>%
  filter(n_htn >=20)

urban_htn$d_id <- as.factor(urban_htn$d_id)

for (i in levels(urban_htn$d_id)){
  print(paste( i))
  df1 = subset(urban_htn, d_id == i)
  model <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data= df1 )
  icc<-calc.icc(model)
  print(icc)
} 

#urban/smoker

urban_smok <- urban %>%
  filter(n_smok >=20)

urban_smok$d_id <- as.factor(urban_smok$d_id)

for (i in levels(urban_smok$d_id)){
  print(paste( i))
  df1 = subset(urban_smok, d_id == i)
  model <- lmer(csmoke ~ 1 + (1|hh_id), data= df1 )
  icc<-calc.icc(model)
  print(icc)
} 

#urban/overweight
urban_ovw <- urban %>%
  filter(n_ovw >=20)

urban_ovw$d_id <- as.factor(urban_ovw$d_id)

for (i in levels(urban_ovw$d_id)){
  print(paste( i))
  df1 = subset(urban_ovw, d_id == i)
  model <- lmer(overweight23 ~ 1 + (1|hh_id), data= df1 )
  icc<-calc.icc(model)
  print(icc)
} 

#urban/obese

urban_obese <- urban %>%
  filter(n_obese >=20)

urban_obese$d_id <- as.factor(urban_obese$d_id)

for (i in levels(urban_obese$d_id)){
  print(paste( i))
  df1 = subset(urban_obese, d_id == i)
  model <- lmer(obese27 ~ 1 + (1|hh_id), data= df1 )
  icc<-calc.icc(model)
  print(icc)
} 


#import Excel file with saved ICCs from calculations above

ICC_hh_d_diab <- read.csv("ICC_hh_by district_diab.csv")
ICC_hh_d_htn <- read.csv("ICC_hh_by district_htn.csv")
ICC_hh_d_csmok <- read.csv("ICC_hh_by district_csmok.csv")
ICC_hh_d_ovw <- read.csv("ICC_hh_by district_ovw.csv")
ICC_hh_d_obese <- read.csv("ICC_hh_by district_obese.csv")

ICC_hh_d <- left_join(ICC_hh_d_htn, ICC_hh_d_diab)# start with htn, because it includes all districts
ICC_hh_d <- left_join(ICC_hh_d, ICC_hh_d_csmok)
ICC_hh_d <- left_join(ICC_hh_d, ICC_hh_d_ovw)
ICC_hh_d <- left_join(ICC_hh_d, ICC_hh_d_obese)

write.csv(ICC_hh_d, file="ICC_hh_d_figure2.csv") # this was uploaded
ICC_hh_d$d_id <- ICC_hh_d$district

ICC_hh_d$d_id<-as.factor(ICC_hh_d$d_id)
distmean.dat$d_id <-as.factor(distmean.dat$d_id)

distmean.dat <- left_join(distmean.dat, ICC_hh_d) 


#create figure for diabetes

distwealthfig <- distmean.dat%>% 
  ggplot(aes(y=medianai, x=ICC_diabetes_hh)) +
  geom_jitter(aes(y=medianai, x=ICC_diabetes_hh, color=zone, shape=zone), size=1) +
  theme_classic() +
  geom_smooth(method='glm', se= FALSE, color="black") +
  labs(x = "ICC at household level",
       y = "Household wealth index", 
       fill="") +
   facet_wrap(~urban_lab)+ 
  theme(axis.text=element_text(size=16,family="Times"),
        axis.title=element_text(size=16),
        legend.text=element_text(size=16,family="Times"),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20),family="Times"),
        axis.title.y = element_text(margin = margin(r = 20),family="Times"),
        strip.text.x = element_text(size=18, ,family="Times"),#rural/urban
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
 scale_x_continuous(breaks = c(0,0.25, 0.5, 0.75,1), limits=c(0,1.1)) +
scale_y_continuous(breaks = c((-2), (-1), 0, 1, 2), limits=c((-3), 3)) +
  coord_fixed(1.1/5,expand=T)
distwealthfig
# face="bold"

#create figure for hypertension

distwealthfig_htn <- distmean.dat %>% 
  ggplot(aes(y=medianai, x=ICC_diabetes_hh)) +
  geom_jitter(aes(y=medianai, x=ICC_diabetes_hh, color=zone, shape=zone), size=1) +
  theme_classic() +
  geom_smooth(method='glm', se= FALSE, color="black") +
  labs(x = "ICC at household level",
       y = "Household wealth index", 
       fill="") +
   facet_wrap(~urban_lab)+ 
  theme(axis.text=element_text(size=16,family="Times"),
        axis.title=element_text(size=16),
        legend.text=element_text(size=16,family="Times"),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20),family="Times"),
        axis.title.y = element_text(margin = margin(r = 20),family="Times"),
        strip.text.x = element_text(size=18, ,family="Times"),#rural/urban
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
 scale_x_continuous(breaks = c(0,0.25, 0.5, 0.75,1), limits=c(0,1.1)) +
scale_y_continuous(breaks = c((-2), (-1), 0, 1, 2), limits=c((-3), 3)) +
  coord_fixed(1.1/5,expand=T)
distwealthfig_htn
# face="bold"

#create figure for current smoker

distwealthfig_csmoke <- distmean.dat %>% 
  ggplot(aes(y=medianai, x=ICC_diabetes_hh)) +
  geom_jitter(aes(y=medianai, x=ICC_diabetes_hh, color=zone, shape=zone), size=1) +
  theme_classic() +
  geom_smooth(method='glm', se= FALSE, color="black") +
  labs(x = "ICC at household level",
       y = "Household wealth index",
       fill="") +
   facet_wrap(~urban_lab)+
  theme(axis.text=element_text(size=16,family="Times"),
        axis.title=element_text(size=16),
        legend.text=element_text(size=16,family="Times"),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20),family="Times"),
        axis.title.y = element_text(margin = margin(r = 20),family="Times"),
        strip.text.x = element_text(size=18, ,family="Times"),#rural/urban
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
 scale_x_continuous(breaks = c(0,0.25, 0.5, 0.75,1), limits=c(0,1.1)) +
scale_y_continuous(breaks = c((-2), (-1), 0, 1, 2), limits=c((-3), 3)) +
  coord_fixed(1.1/5,expand=T)
distwealthfig_csmoke
# face="bold"

#create figure for overweight
 
distwealthfig_overweight <- distmean.dat %>% 
 ggplot(aes(y=medianai, x=ICC_diabetes_hh)) +
  geom_jitter(aes(y=medianai, x=ICC_diabetes_hh, color=zone, shape=zone), size=1) +
  theme_classic() +
  geom_smooth(method='glm', se= FALSE, color="black") +
  labs(x = "ICC at household level",
       y = "Household wealth index", 
       fill="") +
   facet_wrap(~urban_lab)+ 
  theme(axis.text=element_text(size=16,family="Times"),
        axis.title=element_text(size=16),
        legend.text=element_text(size=16,family="Times"),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20),family="Times"),
        axis.title.y = element_text(margin = margin(r = 20),family="Times"),
        strip.text.x = element_text(size=18, ,family="Times"),#rural/urban
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
 scale_x_continuous(breaks = c(0,0.25, 0.5, 0.75,1), limits=c(0,1.1)) +
scale_y_continuous(breaks = c((-2), (-1), 0, 1, 2), limits=c((-3), 3)) +
  coord_fixed(1.1/5,expand=T)
distwealthfig_overweight
# face="bold"

#create figure for obesity

distwealthfig_obese <- distmean.dat %>% 
 ggplot(aes(y=medianai, x=ICC_diabetes_hh)) +
  geom_jitter(aes(y=medianai, x=ICC_diabetes_hh, color=zone, shape=zone), size=1) +
  theme_classic() +
  geom_smooth(method='glm', se= FALSE, color="black") +
  labs(x = "ICC at household level",
       y = "Household wealth index", 
       fill="") +
   facet_wrap(~urban_lab)+
  theme(axis.text=element_text(size=16,family="Times"),
        axis.title=element_text(size=16),
        legend.text=element_text(size=16,family="Times"),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20),family="Times"),
        axis.title.y = element_text(margin = margin(r = 20),family="Times"),
        strip.text.x = element_text(size=18, ,family="Times"),#rural/urban
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
 scale_x_continuous(breaks = c(0,0.25, 0.5, 0.75,1), limits=c(0,1.1)) +
scale_y_continuous(breaks = c((-2), (-1), 0, 1, 2), limits=c((-3), 3)) +
  coord_fixed(1.1/5,expand=T)
distwealthfig_obese
# face="bold"

#calculate p-values

distmean.dat1 <- distmean.dat %>%
  filter(urban_lab=="Rural")

distmean.dat2 <- distmean.dat %>%
  filter(urban_lab=="Urban")

#rural
summary(lm(medianai ~ ICC_diabetes_hh, data=distmean.dat1)) 
summary(lm(medianai ~ ICC_htn_hh, data=distmean.dat1))
summary(lm(medianai ~ ICC_csmoke_hh, data=distmean.dat1)) 
summary(lm(medianai ~ ICC_overweight_hh, data=distmean.dat1))
summary(lm(medianai ~ ICC_obese_hh, data=distmean.dat1)) 


#urban
summary(lm(medianai ~ ICC_diabetes_hh, data=distmean.dat2))
summary(lm(medianai ~ ICC_htn_hh, data=distmean.dat2))
summary(lm(medianai ~ ICC_csmoke_hh, data=distmean.dat2))
summary(lm(medianai ~ ICC_overweight_hh, data=distmean.dat2))
summary(lm(medianai ~ ICC_obese_hh, data=distmean.dat2))

```

```{r APPENDIX: Stratifying by Wealth Quintile}
#made sure each household member has the same household wealth quintile 

#divide dataset into 10 datatsets (rural-urban, each wealth quintile)
rural1 <- india %>%
  filter(urban==0, wealth_quintile_rurb==1)

rural2 <- india %>%
  filter(urban==0, wealth_quintile_rurb==2)

rural3 <- india %>%
  filter(urban==0, wealth_quintile_rurb==3)

rural4 <- india %>%
  filter(urban==0, wealth_quintile_rurb==4)

rural5 <- india %>%
  filter(urban==0, wealth_quintile_rurb==5)

urban1 <- india %>%
  filter(urban==1, wealth_quintile_rurb==1)

urban2 <- india %>%
  filter(urban==1, wealth_quintile_rurb==2)

urban3 <- india %>%
  filter(urban==1, wealth_quintile_rurb==3)

urban4 <- india %>%
  filter(urban==1, wealth_quintile_rurb==4)

urban5 <- india %>%
  filter(urban==1, wealth_quintile_rurb==5)
```

```{r Clustering of outcomes in households}
#take ten divided datasets (rural-urban, wealth quintile 1-5) and calculate for each ICC at household level 
#fixed effect number of people in a household:adultshh has 736777 missings(39%), so I do not use it as a fixed effect

diab_hh_r1 <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=rural1)
calc.icc(diab_hh_r1)
diab_hh_r2 <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=rural2)
calc.icc(diab_hh_r2)
diab_hh_r3 <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=rural3)
calc.icc(diab_hh_r3)
diab_hh_r4 <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=rural4)
calc.icc(diab_hh_r4)
diab_hh_r5 <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=rural5)
calc.icc(diab_hh_r5)
diab_hh_u1 <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=urban1)
calc.icc(diab_hh_u1)
diab_hh_u2 <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=urban2)
calc.icc(diab_hh_u2)
diab_hh_u3 <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=urban3)
calc.icc(diab_hh_u3)
diab_hh_u4 <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=urban4)
calc.icc(diab_hh_u4)
diab_hh_u5 <- lmer(ex_diab_narrow_ind ~ 1 + (1|hh_id), data=urban5)
calc.icc(diab_hh_u5)

htn_hh_r1 <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=rural1)
calc.icc(htn_hh_r1)
htn_hh_r2 <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=rural2)
calc.icc(htn_hh_r2)
htn_hh_r3 <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=rural3)
calc.icc(htn_hh_r3)
htn_hh_r4 <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=rural4)
calc.icc(htn_hh_r4)
htn_hh_r5 <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=rural5)
calc.icc(htn_hh_r5)
htn_hh_u1 <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=urban1)
calc.icc(htn_hh_u1)
htn_hh_u2 <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=urban2)
calc.icc(htn_hh_u2)
htn_hh_u3 <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=urban3)
calc.icc(htn_hh_u3)
htn_hh_u4 <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=urban4)
calc.icc(htn_hh_u4)
htn_hh_u5 <- lmer(ex_htn_narrow_ind ~ 1 + (1|hh_id), data=urban5)
calc.icc(htn_hh_u5)

csmoke_hh_r1 <- lmer(csmoke ~ 1 + (1|hh_id), data=rural1)
calc.icc(csmoke_hh_r1)
csmoke_hh_r2 <- lmer(csmoke ~ 1 + (1|hh_id), data=rural2)
calc.icc(csmoke_hh_r2)
csmoke_hh_r3 <- lmer(csmoke ~ 1 + (1|hh_id), data=rural3)
calc.icc(csmoke_hh_r3)
csmoke_hh_r4 <- lmer(csmoke ~ 1 + (1|hh_id), data=rural4)
calc.icc(csmoke_hh_r4)
csmoke_hh_r5 <- lmer(csmoke ~ 1 + (1|hh_id), data=rural5)
calc.icc(csmoke_hh_r5)
csmoke_hh_u1 <- lmer(csmoke ~ 1 + (1|hh_id), data=urban1)
calc.icc(csmoke_hh_u1)
csmoke_hh_u2 <- lmer(csmoke ~ 1 + (1|hh_id), data=urban2)
calc.icc(csmoke_hh_u2)
csmoke_hh_u3 <- lmer(csmoke ~ 1 + (1|hh_id), data=urban3)
calc.icc(csmoke_hh_u3)
csmoke_hh_u4 <- lmer(csmoke ~ 1 + (1|hh_id), data=urban4)
calc.icc(csmoke_hh_u4)
csmoke_hh_u5 <- lmer(csmoke ~ 1 + (1|hh_id), data=urban5)
calc.icc(csmoke_hh_u5)

overweight_hh_r1 <- lmer(overweight23 ~ 1 + (1|hh_id), data=rural1)
calc.icc(overweight_hh_r1)
overweight_hh_r2 <- lmer(overweight23 ~ 1 + (1|hh_id), data=rural2)
calc.icc(overweight_hh_r2)
overweight_hh_r3 <- lmer(overweight23 ~ 1 + (1|hh_id), data=rural3)
calc.icc(overweight_hh_r3)
overweight_hh_r4 <- lmer(overweight23 ~ 1 + (1|hh_id), data=rural4)
calc.icc(overweight_hh_r4)
overweight_hh_r5 <- lmer(overweight23 ~ 1 + (1|hh_id), data=rural5)
calc.icc(overweight_hh_r5)
overweight_hh_u1 <- lmer(overweight23 ~ 1 + (1|hh_id), data=urban1)
calc.icc(overweight_hh_u1)
overweight_hh_u2 <- lmer(overweight23 ~ 1 + (1|hh_id), data=urban2)
calc.icc(overweight_hh_u2)
overweight_hh_u3 <- lmer(overweight23 ~ 1 + (1|hh_id), data=urban3)
calc.icc(overweight_hh_u3)
overweight_hh_u4 <- lmer(overweight23 ~ 1 + (1|hh_id), data=urban4)
calc.icc(overweight_hh_u4)
overweight_hh_u5 <- lmer(overweight23 ~ 1 + (1|hh_id), data=urban5)
calc.icc(overweight_hh_u5)

obese27_hh_r1 <- lmer(obese27 ~ 1 + (1|hh_id), data=rural1)
calc.icc(obese27_hh_r1)
obese27_hh_r2 <-lmer(obese27 ~ 1 + (1|hh_id), data=rural2)
calc.icc(obese27_hh_r2)
obese27_hh_r3 <- lmer(obese27 ~ 1 + (1|hh_id), data=rural3)
calc.icc(obese27_hh_r3)
obese27_hh_r4 <- lmer(obese27 ~ 1 + (1|hh_id), data=rural4)
calc.icc(obese27_hh_r4)
obese27_hh_r5 <- lmer(obese27 ~ 1 + (1|hh_id), data=rural5)
calc.icc(obese27_hh_r5)
obese27_hh_u1 <- lmer(obese27 ~ 1 + (1|hh_id), data=urban1)
calc.icc(obese27_hh_u1)
obese27_hh_u2 <- lmer(obese27 ~ 1 + (1|hh_id), data=urban2)
calc.icc(obese27_hh_u2)
obese27_hh_u3 <- lmer(obese27 ~ 1 + (1|hh_id), data=urban3)
calc.icc(obese27_hh_u3)
obese27_hh_u4 <- lmer(obese27 ~ 1 + (1|hh_id), data=urban4)
calc.icc(obese27_hh_u4)
obese27_hh_u5 <- lmer(obese27 ~ 1 + (1|hh_id), data=urban5)
calc.icc(obese27_hh_u5)

#calculcate confidence intervals via bootstrapping
#Calculate the bootstrap distribution
boot_diab_hh_r1 <- bootMer(diab_hh_r1, calc.icc, nsim=500)
boot_diab_hh_r2 <- bootMer(diab_hh_r2, calc.icc, nsim=500) 
boot_diab_hh_r3 <- bootMer(diab_hh_r3, calc.icc, nsim=500) 
boot_diab_hh_r4 <- bootMer(diab_hh_r4, calc.icc, nsim=500)
boot_diab_hh_r5 <- bootMer(diab_hh_r5, calc.icc, nsim=500)

boot_diab_hh_u1 <- bootMer(diab_hh_u1, calc.icc, nsim=500)
boot_diab_hh_u2 <- bootMer(diab_hh_u2, calc.icc, nsim=500)
boot_diab_hh_u3 <- bootMer(diab_hh_u3, calc.icc, nsim=500)
boot_diab_hh_u4 <- bootMer(diab_hh_u4, calc.icc, nsim=500)
boot_diab_hh_u5 <- bootMer(diab_hh_u5, calc.icc, nsim=500)

boot_htn_hh_r1 <- bootMer(htn_hh_r1, calc.icc, nsim=500) 
boot_htn_hh_r2 <- bootMer(htn_hh_r2, calc.icc, nsim=500) 
boot_htn_hh_r3 <- bootMer(htn_hh_r3, calc.icc, nsim=500) 
boot_htn_hh_r4 <- bootMer(htn_hh_r4, calc.icc, nsim=500) 
boot_htn_hh_r5 <- bootMer(htn_hh_r5, calc.icc, nsim=500) 
boot_htn_hh_u1 <- bootMer(htn_hh_u1, calc.icc, nsim=500)
boot_htn_hh_u2 <- bootMer(htn_hh_u2, calc.icc, nsim=500)
boot_htn_hh_u3 <- bootMer(htn_hh_u3, calc.icc, nsim=500)
boot_htn_hh_u4 <- bootMer(htn_hh_u4, calc.icc, nsim=500)
boot_htn_hh_u5 <- bootMer(htn_hh_u5, calc.icc, nsim=500)

boot_csmoke_hh_r1 <- bootMer(csmoke_hh_r1, calc.icc, nsim=500) 
boot_csmoke_hh_r2 <- bootMer(csmoke_hh_r2, calc.icc, nsim=500) 
boot_csmoke_hh_r3 <- bootMer(csmoke_hh_r3, calc.icc, nsim=500) 
boot_csmoke_hh_r4 <- bootMer(csmoke_hh_r4, calc.icc, nsim=500) 
boot_csmoke_hh_r5 <- bootMer(csmoke_hh_r5, calc.icc, nsim=500) 
boot_csmoke_hh_u1 <- bootMer(csmoke_hh_u1, calc.icc, nsim=500)
boot_csmoke_hh_u2 <- bootMer(csmoke_hh_u2, calc.icc, nsim=500)
boot_csmoke_hh_u3 <- bootMer(csmoke_hh_u3, calc.icc, nsim=500)
boot_csmoke_hh_u4 <- bootMer(csmoke_hh_u4, calc.icc, nsim=500)
boot_csmoke_hh_u5 <- bootMer(csmoke_hh_u5, calc.icc, nsim=500)

boot_obese27_hh_r1 <- bootMer(obese27_hh_r1, calc.icc, nsim=500) 
boot_obese27_hh_r2 <- bootMer(obese27_hh_r2, calc.icc, nsim=500) 
boot_obese27_hh_r3 <- bootMer(obese27_hh_r3, calc.icc, nsim=500) 
boot_obese27_hh_r4 <- bootMer(obese27_hh_r4, calc.icc, nsim=500) 
boot_obese27_hh_r5 <- bootMer(obese27_hh_r5, calc.icc, nsim=500) 
boot_obese27_hh_u1 <- bootMer(obese27_hh_u1, calc.icc, nsim=500)
boot_obese27_hh_u2 <- bootMer(obese27_hh_u2, calc.icc, nsim=500)
boot_obese27_hh_u3 <- bootMer(obese27_hh_u3, calc.icc, nsim=500)
boot_obese27_hh_u4 <- bootMer(obese27_hh_u4, calc.icc, nsim=500)
boot_obese27_hh_u5 <- bootMer(obese27_hh_u5, calc.icc, nsim=500)

boot_overweight_hh_r1 <- bootMer(overweight_hh_r1, calc.icc, nsim=500) 
boot_overweight_hh_r2 <- bootMer(overweight_hh_r2, calc.icc, nsim=500) 
boot_overweight_hh_r3 <- bootMer(overweight_hh_r3, calc.icc, nsim=500) 
boot_overweight_hh_r4 <- bootMer(overweight_hh_r4, calc.icc, nsim=500) 
boot_overweight_hh_r5 <- bootMer(overweight_hh_r5, calc.icc, nsim=500) 
boot_overweight_hh_u1 <- bootMer(overweight_hh_u1, calc.icc, nsim=500)
boot_overweight_hh_u2 <- bootMer(overweight_hh_u2, calc.icc, nsim=500)
boot_overweight_hh_u3 <- bootMer(overweight_hh_u3, calc.icc, nsim=500)
boot_overweight_hh_u4 <- bootMer(overweight_hh_u4, calc.icc, nsim=500)
boot_overweight_hh_u5 <- bootMer(overweight_hh_u5, calc.icc, nsim=500)


#Draw from the bootstrap distribution the usual 95% upper and lower confidence limits
quantile(boot_diab_hh_r1$t, c(0.025, 0.975))
quantile(boot_diab_hh_r2$t, c(0.025, 0.975))
quantile(boot_diab_hh_r3$t, c(0.025, 0.975))
quantile(boot_diab_hh_r4$t, c(0.025, 0.975))
quantile(boot_diab_hh_r5$t, c(0.025, 0.975))
quantile(boot_diab_hh_u1$t, c(0.025, 0.975))
quantile(boot_diab_hh_u2$t, c(0.025, 0.975))
quantile(boot_diab_hh_u3$t, c(0.025, 0.975))
quantile(boot_diab_hh_u4$t, c(0.025, 0.975))
quantile(boot_diab_hh_u5$t, c(0.025, 0.975))

quantile(boot_htn_hh_r1$t, c(0.025, 0.975))
quantile(boot_htn_hh_r2$t, c(0.025, 0.975))
quantile(boot_htn_hh_r3$t, c(0.025, 0.975))
quantile(boot_htn_hh_r4$t, c(0.025, 0.975))
quantile(boot_htn_hh_r5$t, c(0.025, 0.975))
quantile(boot_htn_hh_u1$t, c(0.025, 0.975))
quantile(boot_htn_hh_u2$t, c(0.025, 0.975))
quantile(boot_htn_hh_u3$t, c(0.025, 0.975))
quantile(boot_htn_hh_u4$t, c(0.025, 0.975))
quantile(boot_htn_hh_u5$t, c(0.025, 0.975))

quantile(boot_csmoke_hh_r1$t, c(0.025, 0.975))
quantile(boot_csmoke_hh_r2$t, c(0.025, 0.975))
quantile(boot_csmoke_hh_r3$t, c(0.025, 0.975))
quantile(boot_csmoke_hh_r4$t, c(0.025, 0.975))
quantile(boot_csmoke_hh_r5$t, c(0.025, 0.975))
quantile(boot_csmoke_hh_u1$t, c(0.025, 0.975))
quantile(boot_csmoke_hh_u2$t, c(0.025, 0.975))
quantile(boot_csmoke_hh_u3$t, c(0.025, 0.975))
quantile(boot_csmoke_hh_u4$t, c(0.025, 0.975))
quantile(boot_csmoke_hh_u5$t, c(0.025, 0.975))

quantile(boot_overweight_hh_r1$t, c(0.025, 0.975))
quantile(boot_overweight_hh_r2$t, c(0.025, 0.975))
quantile(boot_overweight_hh_r3$t, c(0.025, 0.975))
quantile(boot_overweight_hh_r4$t, c(0.025, 0.975))
quantile(boot_overweight_hh_r5$t, c(0.025, 0.975))
quantile(boot_overweight_hh_u1$t, c(0.025, 0.975))
quantile(boot_overweight_hh_u2$t, c(0.025, 0.975))
quantile(boot_overweight_hh_u3$t, c(0.025, 0.975))
quantile(boot_overweight_hh_u4$t, c(0.025, 0.975))
quantile(boot_overweight_hh_u5$t, c(0.025, 0.975))

quantile(boot_obese27_hh_r1$t, c(0.025, 0.975))
quantile(boot_obese27_hh_r2$t, c(0.025, 0.975))
quantile(boot_obese27_hh_r3$t, c(0.025, 0.975))
quantile(boot_obese27_hh_r4$t, c(0.025, 0.975))
quantile(boot_obese27_hh_r5$t, c(0.025, 0.975))
quantile(boot_obese27_hh_u1$t, c(0.025, 0.975))
quantile(boot_obese27_hh_u2$t, c(0.025, 0.975))
quantile(boot_obese27_hh_u3$t, c(0.025, 0.975))
quantile(boot_obese27_hh_u4$t, c(0.025, 0.975))
quantile(boot_obese27_hh_u5$t, c(0.025, 0.975))

```

```{r Clustering in  communities}
#divide dataset into rural-urban, then calculate median of asset index in each psu and divide into quintiles
rural <- india %>%
  filter(urban==0)

urban <-india %>%
  filter(urban==1)

rural <- rural %>% 
  group_by(psuid) %>% 
  mutate(medianai = median(ai_NAT_rurb, na.rm=TRUE)) %>% 
  ungroup()

urban<- urban %>% 
  group_by(psuid) %>% 
  mutate(medianai = median(ai_NAT_rurb, na.rm=TRUE)) %>% 
  ungroup()

# Divide medianai into quintiles
rural <- rural %>% 
  mutate(q_medianai = as.factor(ntile(medianai, 5)))

urban <- urban%>% 
  mutate(q_medianai = as.factor(ntile(medianai, 5)))

rural$wealth_quintile_PSU <-rural$q_medianai
urban$wealth_quintile_PSU <-urban$q_medianai

#divide dataset into 10 datasets by psu wealth quintile
rural1psu <- rural %>%
  filter(wealth_quintile_PSU==1)

rural2psu <- rural %>%
  filter(wealth_quintile_PSU==2)

rural3psu <- rural %>%
  filter(wealth_quintile_PSU==3)

rural4psu <- rural %>%
  filter(wealth_quintile_PSU==4)

rural5psu <- rural %>%
  filter(wealth_quintile_PSU==5)

urban1psu <- urban %>%
  filter(wealth_quintile_PSU==1)

urban2psu <- urban %>%
  filter(wealth_quintile_PSU==2)

urban3psu <- urban %>%
  filter(wealth_quintile_PSU==3)

urban4psu <- urban %>%
  filter(wealth_quintile_PSU==4)

urban5psu <- urban %>%
  filter(wealth_quintile_PSU==5)


###########district level################
#convert d_id
#india$d_id <- as.factor(india$d_id)

rural <- rural %>% 
  group_by(d_id) %>% 
  mutate(medianai_d = median(ai_NAT_rurb, na.rm=TRUE)) %>% 
  ungroup()

urban<- urban %>% 
  group_by(d_id) %>% 
  mutate(medianai_d = median(ai_NAT_rurb, na.rm=TRUE)) %>% 
  ungroup()

# Divide medianai into quintiles
rural <- rural %>% 
  mutate(q_medianai_d = as.factor(ntile(medianai_d, 5)))

urban <- urban%>% 
  mutate(q_medianai_d = as.factor(ntile(medianai_d, 5)))

rural$wealth_quintile_district <-rural$q_medianai_d
urban$wealth_quintile_district <-urban$q_medianai_d

#divide into 10 different datasets by district wealth quintile
rural1district <- rural %>%
  filter(wealth_quintile_district==1)

rural2district <- rural%>%
  filter(wealth_quintile_district==2)

rural3district <- rural %>%
  filter(wealth_quintile_district==3)

rural4district <- rural %>%
  filter(wealth_quintile_district==4)

rural5district <- rural %>%
  filter(wealth_quintile_district==5)

urban1district <- urban %>%
  filter(wealth_quintile_district==1)

urban2district <- urban %>%
  filter( wealth_quintile_district==2)

urban3district <- urban %>%
  filter( wealth_quintile_district==3)

urban4district <- urban %>%
  filter( wealth_quintile_district==4)

urban5district <- urban %>%
  filter(wealth_quintile_district==5)


####################calculate ICCs##################

#on the community level
diab_psu_r1 <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=rural1psu)
calc.icc(diab_psu_r1)
diab_psu_r2 <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=rural2psu)
calc.icc(diab_psu_r2)
diab_psu_r3 <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=rural3psu)
calc.icc(diab_psu_r3)
diab_psu_r4 <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=rural4psu)
calc.icc(diab_psu_r4)
diab_psu_r5 <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=rural5psu)
calc.icc(diab_psu_r5)
diab_psu_u1 <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=urban1psu)
calc.icc(diab_psu_u1)
diab_psu_u2 <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=urban2psu)
calc.icc(diab_psu_u2)
diab_psu_u3 <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=urban3psu)
calc.icc(diab_psu_u3)
diab_psu_u4 <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=urban4psu)
calc.icc(diab_psu_u4)
diab_psu_u5 <- lmer(ex_diab_narrow_ind ~ 1 + (1|psuid), data=urban5psu)
calc.icc(diab_psu_u5)

htn_psu_r1 <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=rural1psu)
calc.icc(htn_psu_r1)
htn_psu_r2 <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=rural2psu)
calc.icc(htn_psu_r2)
htn_psu_r3 <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=rural3psu)
calc.icc(htn_psu_r3)
htn_psu_r4 <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=rural4psu)
calc.icc(htn_psu_r4)
htn_psu_r5 <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=rural5psu)
calc.icc(htn_psu_r5)
htn_psu_u1 <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=urban1psu)
calc.icc(htn_psu_u1)
htn_psu_u2 <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=urban2psu)
calc.icc(htn_psu_u2)
htn_psu_u3 <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=urban3psu)
calc.icc(htn_psu_u3)
htn_psu_u4 <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=urban4psu)
calc.icc(htn_psu_u4)
htn_psu_u5 <- lmer(ex_htn_narrow_ind ~ 1 + (1|psuid), data=urban5psu)
calc.icc(htn_psu_u5)

csmoke_psu_r1 <- lmer(csmoke ~ 1 + (1|psuid), data=rural1psu)
calc.icc(csmoke_psu_r1)
csmoke_psu_r2 <- lmer(csmoke ~ 1 + (1|psuid), data=rural2psu)
calc.icc(csmoke_psu_r2)
csmoke_psu_r3 <- lmer(csmoke ~ 1 + (1|psuid), data=rural3psu)
calc.icc(csmoke_psu_r3)
csmoke_psu_r4 <- lmer(csmoke ~ 1 + (1|psuid), data=rural4psu)
calc.icc(csmoke_psu_r4)
csmoke_psu_r5 <- lmer(csmoke ~ 1 + (1|psuid), data=rural5psu)
calc.icc(csmoke_psu_r5)
csmoke_psu_u1 <- lmer(csmoke ~ 1 + (1|psuid), data=urban1psu)
calc.icc(csmoke_psu_u1)
csmoke_psu_u2 <- lmer(csmoke ~ 1 + (1|psuid), data=urban2psu)
calc.icc(csmoke_psu_u2)
csmoke_psu_u3 <- lmer(csmoke ~ 1 + (1|psuid), data=urban3psu)
calc.icc(csmoke_psu_u3)
csmoke_psu_u4 <- lmer(csmoke ~ 1 + (1|psuid), data=urban4psu)
calc.icc(csmoke_psu_u4)
csmoke_psu_u5 <- lmer(csmoke ~ 1 + (1|psuid), data=urban5psu)
calc.icc(csmoke_psu_u5)

overweight_psu_r1 <- lmer(overweight23 ~ 1 + (1|psuid), data=rural1psu)
calc.icc(overweight_psu_r1)
overweight_psu_r2 <- lmer(overweight23 ~ 1 + (1|psuid), data=rural2psu)
calc.icc(overweight_psu_r2)
overweight_psu_r3 <- lmer(overweight23 ~ 1 + (1|psuid), data=rural3psu)
calc.icc(overweight_psu_r3)
overweight_psu_r4 <- lmer(overweight23 ~ 1 + (1|psuid), data=rural4psu)
calc.icc(overweight_psu_r4)
overweight_psu_r5 <- lmer(overweight23 ~ 1 + (1|psuid), data=rural5psu)
calc.icc(overweight_psu_r5)
overweight_psu_u1 <- lmer(overweight23 ~ 1 + (1|psuid), data=urban1psu)
calc.icc(overweight_psu_u1)
overweight_psu_u2 <- lmer(overweight23 ~ 1 + (1|psuid), data=urban2psu)
calc.icc(overweight_psu_u2)
overweight_psu_u3 <- lmer(overweight23 ~ 1 + (1|psuid), data=urban3psu)
calc.icc(overweight_psu_u3)
overweight_psu_u4 <- lmer(overweight23 ~ 1 + (1|psuid), data=urban4psu)
calc.icc(overweight_psu_u4)
overweight_psu_u5 <- lmer(overweight23 ~ 1 + (1|psuid), data=urban5psu)
calc.icc(overweight_psu_u5)

obese27_psu_r1 <- lmer(obese27 ~ 1 + (1|psuid), data=rural1psu)
calc.icc(obese27_psu_r1)
obese27_psu_r2 <- lmer(obese27 ~ 1 + (1|psuid), data=rural2psu)
calc.icc(obese27_psu_r2)
obese27_psu_r3 <- lmer(obese27 ~ 1 + (1|psuid), data=rural3psu)
calc.icc(obese27_psu_r3)
obese27_psu_r4 <- lmer(obese27 ~ 1 + (1|psuid), data=rural4psu)
calc.icc(obese27_psu_r4)
obese27_psu_r5 <- lmer(obese27 ~ 1 + (1|psuid), data=rural5psu)
calc.icc(obese27_psu_r5)
obese27_psu_u1 <- lmer(obese27 ~ 1 + (1|psuid), data=urban1psu)
calc.icc(obese27_psu_u1)
obese27_psu_u2 <- lmer(obese27 ~ 1 + (1|psuid), data=urban2psu)
calc.icc(obese27_psu_u2)
obese27_psu_u3 <- lmer(obese27 ~ 1 + (1|psuid), data=urban3psu)
calc.icc(obese27_psu_u3)
obese27_psu_u4 <- lmer(obese27 ~ 1 + (1|psuid), data=urban4psu)
calc.icc(obese27_psu_u4)
obese27_psu_u5 <- lmer(obese27 ~ 1 + (1|psuid), data=urban5psu)
calc.icc(obese27_psu_u5)


#calculcate confidence intervals via bootstrapping
#Calculate the bootstrap distribution
boot_diab_psu_r1 <- bootMer(diab_psu_r1, calc.icc, nsim=500) 
boot_diab_psu_r2 <- bootMer(diab_psu_r2, calc.icc, nsim=500) 
boot_diab_psu_r3 <- bootMer(diab_psu_r3, calc.icc, nsim=500) 
boot_diab_psu_r4 <- bootMer(diab_psu_r4, calc.icc, nsim=500) 
boot_diab_psu_r5 <- bootMer(diab_psu_r5, calc.icc, nsim=500) 
boot_diab_psu_u1 <- bootMer(diab_psu_u1, calc.icc, nsim=500)
boot_diab_psu_u2 <- bootMer(diab_psu_u2, calc.icc, nsim=500)
boot_diab_psu_u3 <- bootMer(diab_psu_u3, calc.icc, nsim=500)
boot_diab_psu_u4 <- bootMer(diab_psu_u4, calc.icc, nsim=500)
boot_diab_psu_u5 <- bootMer(diab_psu_u5, calc.icc, nsim=500)

boot_htn_psu_r1 <- bootMer(htn_psu_r1, calc.icc, nsim=500) 
boot_htn_psu_r2 <- bootMer(htn_psu_r2, calc.icc, nsim=500) 
boot_htn_psu_r3 <- bootMer(htn_psu_r3, calc.icc, nsim=500) 
boot_htn_psu_r4 <- bootMer(htn_psu_r4, calc.icc, nsim=500) 
boot_htn_psu_r5 <- bootMer(htn_psu_r5, calc.icc, nsim=500) 
boot_htn_psu_u1 <- bootMer(htn_psu_u1, calc.icc, nsim=500)
boot_htn_psu_u2 <- bootMer(htn_psu_u2, calc.icc, nsim=500)
boot_htn_psu_u3 <- bootMer(htn_psu_u3, calc.icc, nsim=500)
boot_htn_psu_u4 <- bootMer(htn_psu_u4, calc.icc, nsim=500)
boot_htn_psu_u5 <- bootMer(htn_psu_u5, calc.icc, nsim=500)

boot_csmoke_psu_r1 <- bootMer(csmoke_psu_r1, calc.icc, nsim=500) 
boot_csmoke_psu_r2 <- bootMer(csmoke_psu_r2, calc.icc, nsim=500) 
boot_csmoke_psu_r3 <- bootMer(csmoke_psu_r3, calc.icc, nsim=500) 
boot_csmoke_psu_r4 <- bootMer(csmoke_psu_r4, calc.icc, nsim=500) 
boot_csmoke_psu_r5 <- bootMer(csmoke_psu_r5, calc.icc, nsim=500) 
boot_csmoke_psu_u1 <- bootMer(csmoke_psu_u1, calc.icc, nsim=500)
boot_csmoke_psu_u2 <- bootMer(csmoke_psu_u2, calc.icc, nsim=500)
boot_csmoke_psu_u3 <- bootMer(csmoke_psu_u3, calc.icc, nsim=500)
boot_csmoke_psu_u4 <- bootMer(csmoke_psu_u4, calc.icc, nsim=500)
boot_csmoke_psu_u5 <- bootMer(csmoke_psu_u5, calc.icc, nsim=500)

boot_obese27_psu_r1 <- bootMer(obese27_psu_r1, calc.icc, nsim=500) 
boot_obese27_psu_r2 <- bootMer(obese27_psu_r2, calc.icc, nsim=500) 
boot_obese27_psu_r3 <- bootMer(obese27_psu_r3, calc.icc, nsim=500) 
boot_obese27_psu_r4 <- bootMer(obese27_psu_r4, calc.icc, nsim=500) 
boot_obese27_psu_r5 <- bootMer(obese27_psu_r5, calc.icc, nsim=500) 
boot_obese27_psu_u1 <- bootMer(obese27_psu_u1, calc.icc, nsim=500)
boot_obese27_psu_u2 <- bootMer(obese27_psu_u2, calc.icc, nsim=500)
boot_obese27_psu_u3 <- bootMer(obese27_psu_u3, calc.icc, nsim=500)
boot_obese27_psu_u4 <- bootMer(obese27_psu_u4, calc.icc, nsim=500)
boot_obese27_psu_u5 <- bootMer(obese27_psu_u5, calc.icc, nsim=500)

boot_overweight_psu_r1 <- bootMer(overweight_psu_r1, calc.icc, nsim=500) 
boot_overweight_psu_r2 <- bootMer(overweight_psu_r2, calc.icc, nsim=500) 
boot_overweight_psu_r3 <- bootMer(overweight_psu_r3, calc.icc, nsim=500) 
boot_overweight_psu_r4 <- bootMer(overweight_psu_r4, calc.icc, nsim=500) 
boot_overweight_psu_r5 <- bootMer(overweight_psu_r5, calc.icc, nsim=500) 
boot_overweight_psu_u1 <- bootMer(overweight_psu_u1, calc.icc, nsim=500)
boot_overweight_psu_u2 <- bootMer(overweight_psu_u2, calc.icc, nsim=500)
boot_overweight_psu_u3 <- bootMer(overweight_psu_u3, calc.icc, nsim=500)
boot_overweight_psu_u4 <- bootMer(overweight_psu_u4, calc.icc, nsim=500)
boot_overweight_psu_u5 <- bootMer(overweight_psu_u5, calc.icc, nsim=500)


#Draw from the bootstrap distribution the usual 95% upper and lower confidence limits
quantile(boot_diab_psu_r1$t, c(0.025, 0.975))
quantile(boot_diab_psu_r2$t, c(0.025, 0.975))
quantile(boot_diab_psu_r3$t, c(0.025, 0.975))
quantile(boot_diab_psu_r4$t, c(0.025, 0.975))
quantile(boot_diab_psu_r5$t, c(0.025, 0.975))
quantile(boot_diab_psu_u1$t, c(0.025, 0.975))
quantile(boot_diab_psu_u2$t, c(0.025, 0.975))
quantile(boot_diab_psu_u3$t, c(0.025, 0.975))
quantile(boot_diab_psu_u4$t, c(0.025, 0.975))
quantile(boot_diab_psu_u5$t, c(0.025, 0.975))

quantile(boot_htn_psu_r1$t, c(0.025, 0.975))
quantile(boot_htn_psu_r2$t, c(0.025, 0.975))
quantile(boot_htn_psu_r3$t, c(0.025, 0.975))
quantile(boot_htn_psu_r4$t, c(0.025, 0.975))
quantile(boot_htn_psu_r5$t, c(0.025, 0.975))
quantile(boot_htn_psu_u1$t, c(0.025, 0.975))
quantile(boot_htn_psu_u2$t, c(0.025, 0.975))
quantile(boot_htn_psu_u3$t, c(0.025, 0.975))
quantile(boot_htn_psu_u4$t, c(0.025, 0.975))
quantile(boot_htn_psu_u5$t, c(0.025, 0.975))

quantile(boot_csmoke_psu_r1$t, c(0.025, 0.975))
quantile(boot_csmoke_psu_r2$t, c(0.025, 0.975))
quantile(boot_csmoke_psu_r3$t, c(0.025, 0.975))
quantile(boot_csmoke_psu_r4$t, c(0.025, 0.975))
quantile(boot_csmoke_psu_r5$t, c(0.025, 0.975))
quantile(boot_csmoke_psu_u1$t, c(0.025, 0.975))
quantile(boot_csmoke_psu_u2$t, c(0.025, 0.975))
quantile(boot_csmoke_psu_u3$t, c(0.025, 0.975))
quantile(boot_csmoke_psu_u4$t, c(0.025, 0.975))
quantile(boot_csmoke_psu_u5$t, c(0.025, 0.975))



quantile(boot_overweight_psu_r1$t, c(0.025, 0.975))
quantile(boot_overweight_psu_r2$t, c(0.025, 0.975))
quantile(boot_overweight_psu_r3$t, c(0.025, 0.975))
quantile(boot_overweight_psu_r4$t, c(0.025, 0.975))
quantile(boot_overweight_psu_r5$t, c(0.025, 0.975))
quantile(boot_overweight_psu_u1$t, c(0.025, 0.975))
quantile(boot_overweight_psu_u2$t, c(0.025, 0.975))
quantile(boot_overweight_psu_u3$t, c(0.025, 0.975))
quantile(boot_overweight_psu_u4$t, c(0.025, 0.975))
quantile(boot_overweight_psu_u5$t, c(0.025, 0.975))

quantile(boot_obese27_psu_r1$t, c(0.025, 0.975))
quantile(boot_obese27_psu_r2$t, c(0.025, 0.975))
quantile(boot_obese27_psu_r3$t, c(0.025, 0.975))
quantile(boot_obese27_psu_r4$t, c(0.025, 0.975))
quantile(boot_obese27_psu_r5$t, c(0.025, 0.975))
quantile(boot_obese27_psu_u1$t, c(0.025, 0.975))
quantile(boot_obese27_psu_u2$t, c(0.025, 0.975))
quantile(boot_obese27_psu_u3$t, c(0.025, 0.975))
quantile(boot_obese27_psu_u4$t, c(0.025, 0.975))
quantile(boot_obese27_psu_u5$t, c(0.025, 0.975))

```

```{r Clustering on the district level}

diab_district_r1 <- lmer(ex_diab_narrow_ind ~ 1 + (1|d_id),data=rural1district)
calc.icc(diab_district_r1)
diab_district_r2 <- lmer(ex_diab_narrow_ind ~ 1 + (1|d_id),data=rural2district)
calc.icc(diab_district_r2)
diab_district_r3 <- lmer(ex_diab_narrow_ind ~ 1 + (1|d_id),data=rural3district)
calc.icc(diab_district_r3)
diab_district_r4 <- lmer(ex_diab_narrow_ind ~ 1 + (1|d_id),data=rural4district)
calc.icc(diab_district_r4)
diab_district_r5 <- lmer(ex_diab_narrow_ind ~ 1 + (1|d_id),data=rural5district)
calc.icc(diab_district_r5)
diab_district_u1 <- lmer(ex_diab_narrow_ind ~ 1 + (1|d_id),data=urban1district)
calc.icc(diab_district_u1)
diab_district_u2 <- lmer(ex_diab_narrow_ind ~ 1 + (1|d_id),data=urban2district)
calc.icc(diab_district_u2)
diab_district_u3 <- lmer(ex_diab_narrow_ind ~ 1 + (1|d_id),data=urban3district)
calc.icc(diab_district_u3)
diab_district_u4 <- lmer(ex_diab_narrow_ind ~ 1 + (1|d_id),data=urban4district)
calc.icc(diab_district_u4)
diab_district_u5 <- lmer(ex_diab_narrow_ind ~ 1 + (1|d_id),data=urban5district)
calc.icc(diab_district_u5)

htn_district_r1 <- lmer(ex_htn_narrow_ind ~ 1 + (1|d_id), data=rural1district)
calc.icc(htn_district_r1)
htn_district_r2 <- lmer(ex_htn_narrow_ind ~ 1 + (1|d_id), data=rural2district)
calc.icc(htn_district_r2)
htn_district_r3 <- lmer(ex_htn_narrow_ind ~ 1 + (1|d_id), data=rural3district)
calc.icc(htn_district_r3)
htn_district_r4 <- lmer(ex_htn_narrow_ind ~ 1 + (1|d_id), data=rural4district)
calc.icc(htn_district_r4)
htn_district_r5 <- lmer(ex_htn_narrow_ind ~ 1 + (1|d_id), data=rural5district)
calc.icc(htn_district_r5)
htn_district_u1 <- lmer(ex_htn_narrow_ind ~ 1 + (1|d_id), data=urban1district)
calc.icc(htn_district_u1)
htn_district_u2 <- lmer(ex_htn_narrow_ind ~ 1 + (1|d_id), data=urban2district)
calc.icc(htn_district_u2)
htn_district_u3 <- lmer(ex_htn_narrow_ind ~ 1 + (1|d_id), data=urban3district)
calc.icc(htn_district_u3)
htn_district_u4 <- lmer(ex_htn_narrow_ind ~ 1 + (1|d_id), data=urban4district)
calc.icc(htn_district_u4)
htn_district_u5 <- lmer(ex_htn_narrow_ind ~ 1 + (1|d_id), data=urban5district)
calc.icc(htn_district_u5)

csmoke_district_r1 <- lmer(csmoke ~ 1 + (1|d_id), data=rural1district)
calc.icc(csmoke_district_r1)
csmoke_district_r2 <- lmer(csmoke ~ 1 + (1|d_id), data=rural2district)
calc.icc(csmoke_district_r2)
csmoke_district_r3 <- lmer(csmoke ~ 1 + (1|d_id), data=rural3district)
calc.icc(csmoke_district_r3)
csmoke_district_r4 <- lmer(csmoke ~ 1 + (1|d_id), data=rural4district)
calc.icc(csmoke_district_r4)
csmoke_district_r5 <- lmer(csmoke ~ 1 + (1|d_id), data=rural5district)
calc.icc(csmoke_district_r5)
csmoke_district_u1 <- lmer(csmoke ~ 1 + (1|d_id), data=urban1district)
calc.icc(csmoke_district_u1)
csmoke_district_u2 <- lmer(csmoke ~ 1 + (1|d_id), data=urban2district)
calc.icc(csmoke_district_u2)
csmoke_district_u3 <- lmer(csmoke ~ 1 + (1|d_id), data=urban3district)
calc.icc(csmoke_district_u3)
csmoke_district_u4 <- lmer(csmoke ~ 1 + (1|d_id), data=urban4district)
calc.icc(csmoke_district_u4)
csmoke_district_u5 <- lmer(csmoke~ 1 + (1|d_id), data=urban5district)
calc.icc(csmoke_district_u5)

overweight_district_r1 <- lmer(overweight23 ~ 1 + (1|d_id), data=rural1district)
calc.icc(overweight_district_r1)
overweight_district_r2 <- lmer(overweight23 ~ 1 + (1|d_id), data=rural2district)
calc.icc(overweight_district_r2)
overweight_district_r3 <- lmer(overweight23 ~ 1 + (1|d_id), data=rural3district)
calc.icc(overweight_district_r3)
overweight_district_r4 <- lmer(overweight23 ~ 1 + (1|d_id), data=rural4district)
calc.icc(overweight_district_r4)
overweight_district_r5 <- lmer(overweight23 ~ 1 + (1|d_id), data=rural5district)
calc.icc(overweight_district_r5)
overweight_district_u1 <- lmer(overweight23 ~ 1 + (1|d_id), data=urban1district)
calc.icc(overweight_district_u1)
overweight_district_u2 <- lmer(overweight23 ~ 1 + (1|d_id), data=urban2district)
calc.icc(overweight_district_u2)
overweight_district_u3 <- lmer(overweight23 ~ 1 + (1|d_id), data=urban3district)
calc.icc(overweight_district_u3)
overweight_district_u4 <- lmer(overweight23 ~ 1 + (1|d_id), data=urban4district)
calc.icc(overweight_district_u4)
overweight_district_u5 <- lmer(overweight23 ~ 1 + (1|d_id), data=urban5district)
calc.icc(overweight_district_u5)

obese27_district_r1 <- lmer(obese27 ~ 1 + (1|d_id), data=rural1district)
calc.icc(obese27_district_r1)
obese27_district_r2 <- lmer(obese27 ~ 1 + (1|d_id), data=rural2district)
calc.icc(obese27_district_r2)
obese27_district_r3 <- lmer(obese27 ~ 1 + (1|d_id), data=rural3district)
calc.icc(obese27_district_r3)
obese27_district_r4 <- lmer(obese27~ 1 + (1|d_id), data=rural4district)
calc.icc(obese27_district_r4)
obese27_district_r5 <- lmer(obese27 ~ 1 + (1|d_id), data=rural5district)
calc.icc(obese27_district_r5)
obese27_district_u1 <- lmer(obese27 ~ 1 + (1|d_id), data=urban1district)
calc.icc(obese27_district_u1)
obese27_district_u2 <- lmer(obese27~ 1 + (1|d_id), data=urban2district)
calc.icc(obese27_district_u2)
obese27_district_u3 <- lmer(obese27 ~ 1 + (1|d_id), data=urban3district)
calc.icc(obese27_district_u3)
obese27_district_u4 <- lmer(obese27 ~ 1 + (1|d_id), data=urban4district)
calc.icc(obese27_district_u4)
obese27_district_u5 <- lmer(obese27 ~ 1 + (1|d_id), data=urban5district)
calc.icc(obese27_district_u5)


#calculate confidence intervals via bootstrapping

boot_diab_district_r1 <- bootMer(diab_district_r1, calc.icc, nsim=500) 
boot_diab_district_r2 <- bootMer(diab_district_r2, calc.icc, nsim=500) 
boot_diab_district_r3 <- bootMer(diab_district_r3, calc.icc, nsim=500) 
boot_diab_district_r4 <- bootMer(diab_district_r4, calc.icc, nsim=500) 
boot_diab_district_r5 <- bootMer(diab_district_r5, calc.icc, nsim=500) 
boot_diab_district_u1 <- bootMer(diab_district_u1, calc.icc, nsim=500)
boot_diab_district_u2 <- bootMer(diab_district_u2, calc.icc, nsim=500)
boot_diab_district_u3 <- bootMer(diab_district_u3, calc.icc, nsim=500)
boot_diab_district_u4 <- bootMer(diab_district_u4, calc.icc, nsim=500)
boot_diab_district_u5 <- bootMer(diab_district_u5, calc.icc, nsim=500)

boot_htn_district_r1 <- bootMer(htn_district_r1, calc.icc, nsim=500) 
boot_htn_district_r2 <- bootMer(htn_district_r2, calc.icc, nsim=500) 
boot_htn_district_r3 <- bootMer(htn_district_r3, calc.icc, nsim=500) 
boot_htn_district_r4 <- bootMer(htn_district_r4, calc.icc, nsim=500) 
boot_htn_district_r5 <- bootMer(htn_district_r5, calc.icc, nsim=500) 
boot_htn_district_u1 <- bootMer(htn_district_u1, calc.icc, nsim=500)
boot_htn_district_u2 <- bootMer(htn_district_u2, calc.icc, nsim=500)
boot_htn_district_u3 <- bootMer(htn_district_u3, calc.icc, nsim=500)
boot_htn_district_u4 <- bootMer(htn_district_u4, calc.icc, nsim=500)
boot_htn_district_u5 <- bootMer(htn_district_u5, calc.icc, nsim=500)

boot_csmoke_district_r1 <- bootMer(csmoke_district_r1, calc.icc, nsim=500) 
boot_csmoke_district_r2 <- bootMer(csmoke_district_r2, calc.icc, nsim=500) 
boot_csmoke_district_r3 <- bootMer(csmoke_district_r3, calc.icc, nsim=500) 
boot_csmoke_district_r4 <- bootMer(csmoke_district_r4, calc.icc, nsim=500) 
boot_csmoke_district_r5 <- bootMer(csmoke_district_r5, calc.icc, nsim=500) 
boot_csmoke_district_u1 <- bootMer(csmoke_district_u1, calc.icc, nsim=500)
boot_csmoke_district_u2 <- bootMer(csmoke_district_u2, calc.icc, nsim=500)
boot_csmoke_district_u3 <- bootMer(csmoke_district_u3, calc.icc, nsim=500)
boot_csmoke_district_u4 <- bootMer(csmoke_district_u4, calc.icc, nsim=500)
boot_csmoke_district_u5 <- bootMer(csmoke_district_u5, calc.icc, nsim=500)

boot_obese27_district_r1 <- bootMer(obese27_district_r1, calc.icc, nsim=500) 
boot_obese27_district_r2 <- bootMer(obese27_district_r2, calc.icc, nsim=500) 
boot_obese27_district_r3 <- bootMer(obese27_district_r3, calc.icc, nsim=500) 
boot_obese27_district_r4 <- bootMer(obese27_district_r4, calc.icc, nsim=500) 
boot_obese27_district_r5 <- bootMer(obese27_district_r5, calc.icc, nsim=500) 
boot_obese27_district_u1 <- bootMer(obese27_district_u1, calc.icc, nsim=500)
boot_obese27_district_u2 <- bootMer(obese27_district_u2, calc.icc, nsim=500)
boot_obese27_district_u3 <- bootMer(obese27_district_u3, calc.icc, nsim=500)
boot_obese27_district_u4 <- bootMer(obese27_district_u4, calc.icc, nsim=500)
boot_obese27_district_u5 <- bootMer(obese27_district_u5, calc.icc, nsim=500)

boot_overweight_district_r1 <- bootMer(overweight_district_r1, calc.icc, nsim=500) 
boot_overweight_district_r2 <- bootMer(overweight_district_r2, calc.icc, nsim=500) 
boot_overweight_district_r3 <- bootMer(overweight_district_r3, calc.icc, nsim=500) 
boot_overweight_district_r4 <- bootMer(overweight_district_r4, calc.icc, nsim=500) 
boot_overweight_district_r5 <- bootMer(overweight_district_r5, calc.icc, nsim=500) 
boot_overweight_district_u1 <- bootMer(overweight_district_u1, calc.icc, nsim=500)
boot_overweight_district_u2 <- bootMer(overweight_district_u2, calc.icc, nsim=500)
boot_overweight_district_u3 <- bootMer(overweight_district_u3, calc.icc, nsim=500)
boot_overweight_district_u4 <- bootMer(overweight_district_u4, calc.icc, nsim=500)
boot_overweight_district_u5 <- bootMer(overweight_district_u5, calc.icc, nsim=500)


#Draw from the bootstrap distribution the usual 95% upper and lower confidence limits
quantile(boot_diab_district_r1$t, c(0.025, 0.975))
quantile(boot_diab_district_r2$t, c(0.025, 0.975))
quantile(boot_diab_district_r3$t, c(0.025, 0.975))
quantile(boot_diab_district_r4$t, c(0.025, 0.975))
quantile(boot_diab_district_r5$t, c(0.025, 0.975))
quantile(boot_diab_district_u1$t, c(0.025, 0.975))
quantile(boot_diab_district_u2$t, c(0.025, 0.975))
quantile(boot_diab_district_u3$t, c(0.025, 0.975))
quantile(boot_diab_district_u4$t, c(0.025, 0.975))
quantile(boot_diab_district_u5$t, c(0.025, 0.975))

quantile(boot_htn_district_r1$t, c(0.025, 0.975))
quantile(boot_htn_district_r2$t, c(0.025, 0.975))
quantile(boot_htn_district_r3$t, c(0.025, 0.975))
quantile(boot_htn_district_r4$t, c(0.025, 0.975))
quantile(boot_htn_district_r5$t, c(0.025, 0.975))
quantile(boot_htn_district_u1$t, c(0.025, 0.975))
quantile(boot_htn_district_u2$t, c(0.025, 0.975))
quantile(boot_htn_district_u3$t, c(0.025, 0.975))
quantile(boot_htn_district_u4$t, c(0.025, 0.975))
quantile(boot_htn_district_u5$t, c(0.025, 0.975))

quantile(boot_csmoke_district_r1$t, c(0.025, 0.975))
quantile(boot_csmoke_district_r2$t, c(0.025, 0.975))
quantile(boot_csmoke_district_r3$t, c(0.025, 0.975))
quantile(boot_csmoke_district_r4$t, c(0.025, 0.975))
quantile(boot_csmoke_district_r5$t, c(0.025, 0.975))
quantile(boot_csmoke_district_u1$t, c(0.025, 0.975))
quantile(boot_csmoke_district_u2$t, c(0.025, 0.975))
quantile(boot_csmoke_district_u3$t, c(0.025, 0.975))
quantile(boot_csmoke_district_u4$t, c(0.025, 0.975))
quantile(boot_csmoke_district_u5$t, c(0.025, 0.975))


quantile(boot_overweight_district_r1$t, c(0.025, 0.975))
quantile(boot_overweight_district_r2$t, c(0.025, 0.975))
quantile(boot_overweight_district_r3$t, c(0.025, 0.975))
quantile(boot_overweight_district_r4$t, c(0.025, 0.975))
quantile(boot_overweight_district_r5$t, c(0.025, 0.975))
quantile(boot_overweight_district_u1$t, c(0.025, 0.975))
quantile(boot_overweight_district_u2$t, c(0.025, 0.975))
quantile(boot_overweight_district_u3$t, c(0.025, 0.975))
quantile(boot_overweight_district_u4$t, c(0.025, 0.975))
quantile(boot_overweight_district_u5$t, c(0.025, 0.975))

quantile(boot_obese27_district_r1$t, c(0.025, 0.975))
quantile(boot_obese27_district_r2$t, c(0.025, 0.975))
quantile(boot_obese27_district_r3$t, c(0.025, 0.975))
quantile(boot_obese27_district_r4$t, c(0.025, 0.975))
quantile(boot_obese27_district_r5$t, c(0.025, 0.975))
quantile(boot_obese27_district_u1$t, c(0.025, 0.975))
quantile(boot_obese27_district_u2$t, c(0.025, 0.975))
quantile(boot_obese27_district_u3$t, c(0.025, 0.975))
quantile(boot_obese27_district_u4$t, c(0.025, 0.975))
quantile(boot_obese27_district_u5$t, c(0.025, 0.975))

```

```{APPENDIX: Figure S1/S2}

#Figure S1 :Regression of CVD risk factor prevalence ~ ICC on household/community level

####### Adding state labels
india <- mutate(india, state_lab = fct_recode(ex_state_ind,
                                                        "HP" = "Himachal Pradesh",
                                                        "PB" = "Punjab",
                                                        "CH" = "Chandigarh",
                                                        "HR" = "Haryana",
                                                        "DL" = "Delhi",
                                                        "SK" = "Sikkim",
                                                        "DD" = "Daman and Diu",
                                                        "AR" = "Arunachal Pradesh",
                                                        "NL" = "Nagaland",
                                                        "MN" = "Manipur",
                                                        "MZ" = "Mizoram",
                                                        "TR" = "Tripura",
                                                        "ML" = "Meghalaya",
                                                        "WB" = "West Bengal",
                                                        "MH" = "Maharashtra",
                                                        "AP" = "Andhra Pradesh",
                                                        "KA" = "Karnataka",
                                                        "GA" = "Goa",
                                                        "KL" = "Kerala",
                                                        "PY" = "Puducherry",
                                                        "TN" = "Tamil Nadu",
                                                        "AN" = "Andaman and Nicobar",
                                                        "TS" = "Telangana",
                                                        "UK" = "Uttarakhand",
                                                        "RJ" = "Rajasthan",
                                                        "UP" = "Uttar Pradesh",
                                                        "BR" = "Bihar",
                                                        "AS" = "Assam",
                                                        "JK" = "Jammu and Kashmir",
                                                        "GJ" = "Gujarat",
                                                        "JH" = "Jharkhand",
                                                        "OD" = "Odisha",
                                                        "CT" = "Chhattisgarh",
                                                        "MP" = "Madhya Pradesh"))



india <- mutate(india,
                     urban_dbl = as.numeric(urban))

#install.packages("spatstat")


library(spatstat)

#calculcate prevalence of CVD risk factors

#india$ex_htn_narrow_ind<-as.numeric(levels(india$ex_htn_narrow_ind))[india_reg$ex_htn_narrow_ind]
#india_reg$ex_diab_narrow_ind<-as.numeric(levels(india$ex_diab_narrow_ind))[india_reg$ex_diab_narrow_ind]
#india$csmoke<-as.numeric(levels(india$csmoke))[india$csmoke]
#india$overweight23<-as.numeric(levels(india$overweight23))[india$overweight23]
#india$obese27<-as.numeric(levels(india$obese27))[india$obese27]

stateprev.dat <- india_reg %>% 
  group_by(ex_state_ind) %>%
  mutate(diab_st=100*weighted.mean(ex_diab_narrow_ind,sworld_weight_india,na.rm=TRUE,vartype="ci")) %>%
  mutate(htn_st=100*weighted.mean(ex_htn_narrow_ind,sworld_weight_india,na.rm=TRUE,vartype="ci")) %>%
  mutate(smok_st=100*weighted.mean(csmoke,sworld_weight_india,na.rm=TRUE,vartype="ci")) %>%
  mutate(obese_st=100*weighted.mean(obese27,sworld_weight_india,na.rm=TRUE,vartype="ci")) %>%
  mutate(overweight_st=100*weighted.mean(overweight23,sworld_weight_india,na.rm=TRUE,vartype="ci")) %>%
  filter(row_number()==1) %>%
  ungroup() %>%
  dplyr::select(state_lab,diab_st,htn_st, smok_st,obese_st,overweight_st,ex_state_ind,zone) 

#add ICCs
ICC_hh <- read.csv("ICC_by_state_dlhsahs_hh.csv") #Excel file
ICC_comm <- read.csv("ICC_by_state_dlhsahs_comm.csv") #Excel file

stateprev.dat <-left_join( stateprev.dat, ICC_hh, by=c("ex_state_ind"="state")) 
stateprev.dat <-left_join( stateprev.dat, ICC_comm, by=c("ex_state_ind"="state")) 


#####Diabetes/household level

statediabhhfig <- stateprev.dat %>%
 
  ggplot(aes(y=diab_st, x=ICC_diabetes_hh)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=diab_st, x=ICC_diabetes_hh, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=diab_st, x=ICC_diabetes_hh, label = state_lab)) +
  theme_classic() +
  labs(x = "ICC at household level",
       y = " Raised BG prevalence, in %",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(10,20,30,40,50), limits=c(0, 60)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/60,expand=F)
statediabhhfig

summary(lm(diab_st ~ ICC_diabetes_hh, data=stateprev.dat))


#######htn/hh

statehtnhhfig <- stateprev.dat %>%
 
  ggplot(aes(y=htn_st, x=ICC_htn_hh)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=htn_st, x=ICC_htn_hh, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=htn_st, x=ICC_htn_hh, label = state_lab)) +
  theme_classic() +
  labs(x = "ICC at household level",
       y = " Raised BP prevalence, in %",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(10,20,30,40,50), limits=c(0, 60)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/60,expand=F)
statehtnhhfig

summary(lm(htn_st ~ ICC_htn_hh, data=stateprev.dat))

###smoke/hh

statesmokehhfig <- stateprev.dat %>%
 
  ggplot(aes(y=smok_st, x=ICC_csmoke_hh)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=smok_st, x=ICC_csmoke_hh, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=smok_st, x=ICC_csmoke_hh, label = state_lab)) +
  theme_classic() +
  labs(x = "ICC at household level",
       y = " Smoking prevalence, in %",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
 scale_y_continuous(breaks = c(0, 10, 20, 30, 40, 50), limits=c(0, 60)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0.000, 0.4)) +
  coord_fixed(0.4/60,expand=F)
statesmokehhfig

summary(lm(smok_st ~ ICC_csmoke_hh, data=stateprev.dat))

##overweight/hh

stateoverweighthhfig <- stateprev.dat %>%
 
  ggplot(aes(y=overweight_st, x=ICC_overweight_hh)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=overweight_st, x=ICC_overweight_hh, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=overweight_st, x=ICC_overweight_hh, label = state_lab)) +
  theme_classic() +
  labs(x = "ICC at household level",
       y = " Overweight prevalence, in %",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(0, 10, 20,30, 40, 50), limits=c(0, 60)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/60,expand=F)
stateoverweighthhfig

summary(lm(overweight_st ~ ICC_overweight_hh, data=stateprev.dat))

##obese/hh

stateobesehhfig <- stateprev.dat %>%
 
  ggplot(aes(y=obese_st, x=ICC_obese_hh)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=obese_st, x=ICC_obese_hh, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=obese_st, x=ICC_obese_hh, label = state_lab)) +
  theme_classic() +
  labs(x = "ICC at household level, in %",
       y = " Obesity prevalence, in %",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(0, 10, 20,30, 40, 50), limits=c(0, 60)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/60,expand=F)
stateobesehhfig

summary(lm(obese_st ~ ICC_obese_hh, data=stateprev.dat))

###########figures for ICC on community level#####

statediabcommfig <- stateprev.dat %>%
 
  ggplot(aes(y=diab_st, x=ICC_diab_comm)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=diab_st, x=ICC_diab_comm, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=diab_st, x=ICC_diab_comm, label = state_lab)) +
  theme_classic() +
  labs(x = "ICC at community level, in %",
       y = " Raised BG prevalence, in %",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(0, 10, 20,30, 40, 50), limits=c(0, 60)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/60,expand=F)
statediabcommfig

#######htn/comm

statehtncommfig <- stateprev.dat %>%
 
  ggplot(aes(y=htn_st, x=ICC_htn_comm)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=htn_st, x=ICC_htn_comm, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=htn_st, x=ICC_htn_comm, label = state_lab)) +
  theme_classic() +
  labs(x = "ICC at community level",
       y = " Raised BP prevalence, in %",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(0,10,20,30,40,50), limits=c(0, 60)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/60,expand=F)
statehtncommfig

###smoke/comm

statesmokecommfig <- stateprev.dat %>%
 
  ggplot(aes(y=smok_st, x=ICC_csmoke_comm)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=smok_st, x=ICC_csmoke_comm, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=smok_st, x=ICC_csmoke_comm, label = state_lab)) +
  theme_classic() +
  labs(x = "ICC at community level, in %",
       y = " Smoking prevalence, in %",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(0,10,20,30,40,50), limits=c(0, 60)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/60,expand=F)
statesmokecommfig

##overweight/comm

stateoverweightcommfig <- stateprev.dat %>%
 
  ggplot(aes(y=overweight_st, x=ICC_overweight_comm)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=overweight_st, x=ICC_overweight_comm, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=overweight_st, x=ICC_overweight_comm, label = state_lab)) +
  theme_classic() +
  labs(x = "ICC at community level, in %",
       y = " Overweight prevalence, in %",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(0,10,20,30,40,50), limits=c(0, 60)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/60,expand=F)
stateoverweightcommfig

##obese/comm

stateobesecommfig <- stateprev.dat %>%
 
  ggplot(aes(y=obese_st, x=ICC_obese_comm)) +
  geom_jitter(aes(y=obese_st, x=ICC_obese_comm, color=zone, shape=zone), size=2.5) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_text_repel(aes(y=obese_st, x=ICC_obese_comm, label = state_lab)) +
  theme_classic() +
  labs(x = "ICC at community level, in %",
       y = " Obesity prevalence, in %",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(0,10,20,30,40,50), limits=c(0, 60)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/60,expand=F)
stateobesecommfig

summary(lm(diab_st ~ ICC_diab_comm, data=stateprev.dat))
summary(lm(htn_st ~ ICC_htn_comm, data=stateprev.dat))
summary(lm(smok_st ~ ICC_csmoke_comm, data=stateprev.dat))
summary(lm(overweight_st ~ ICC_overweight_comm, data=stateprev.dat))
summary(lm(obese_st ~ ICC_obese_comm, data=stateprev.dat))

########Figure S2 :ICC on the household/community level in relation to GDP per capita by state

#import GDP file
#No data is available for the union territories of Dadra and Nagar Haveli, Daman and Diu and Lakshadweep
#State-wise Per Capita Net State Domestic Product (Per Capita NSDP) at CONSTANT PRICES (Rs.) 2013-2014 (Directorate of Economics & Statistics of respective State Governments)				
GDPpercapita<-read.csv("gdppercapita_india.csv")

#add GDP per capita to stateprev.dat
stateprev.dat <-left_join( stateprev.dat,GDPpercapita, by=c("ex_state_ind"="state")) 
stateprev.dat$ex_state_ind <-as.factor(stateprev.dat$ex_state_ind) 

#####Diab/hh
statediabhhfiggdp <- stateprev.dat %>%
 
  ggplot(aes(y=GDP, x=ICC_diabetes_hh)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=GDP, x=ICC_diabetes_hh, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=GDP, x=ICC_diabetes_hh, label = state_lab)) +
  theme_classic() +
  labs(x = "ICC at household level",
       y = " GDP per capita (INR)",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(50000,100000,150000,200000,250000), limits=c(0, 300000)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/300000,expand=F)

statediabhhfiggdp

#calculate p-values
summary(lm(GDP ~ ICC_diabetes_hh, data=stateprev.dat))
summary(lm(GDP ~ ICC_htn_hh, data=stateprev.dat))
summary(lm(GDP ~ ICC_csmoke_hh, data=stateprev.dat))
summary(lm(GDP ~ ICC_overweight_hh, data=stateprev.dat))
summary(lm(GDP ~ ICC_obese_hh, data=stateprev.dat))

#######htn/hh

statehtnhhfiggdp <- stateprev.dat %>%
 
  ggplot(aes(y=GDP, x=ICC_htn_hh)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=GDP, x=ICC_htn_hh, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=GDP, x=ICC_htn_hh, label = state_lab)) +
  theme_classic() +
    labs(x = "ICC at household level",
       y = " GDP per capita (INR)",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(50000,100000,150000,200000,250000), limits=c(0, 300000)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/300000,expand=F)
statehtnhhfiggdp

###smoke/hh

statesmokehhfiggdp <- stateprev.dat %>%
 
  ggplot(aes(y=GDP, x=ICC_csmoke_hh)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=GDP, x=ICC_csmoke_hh, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=GDP, x=ICC_csmoke_hh, label = state_lab)) +
  theme_classic() +
    labs(x = "ICC at household level",
       y = " GDP per capita (INR)",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
 scale_y_continuous(breaks = c(50000,100000,150000,200000,250000), limits=c(0, 300000)) +
 scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed((0.4/300000),expand=T)
statesmokehhfiggdp

##overweight/hh

stateoverweighthhfiggdp <- stateprev.dat %>%
 
  ggplot(aes(y=GDP, x=ICC_overweight_hh)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=GDP, x=ICC_overweight_hh, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=GDP, x=ICC_overweight_hh, label = state_lab)) +
 theme_classic() +
    labs(x = "ICC at household level",
       y = " GDP per capita (INR)",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(50000,100000,150000,200000,250000), limits=c(0, 300000)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/300000,expand=F)
stateoverweighthhfiggdp

##obese/hh

stateobesehhfiggdp <- stateprev.dat %>%
 
  ggplot(aes(y=GDP, x=ICC_obese_hh)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=GDP, x=ICC_obese_hh, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=GDP, x=ICC_obese_hh, label = state_lab)) +
 theme_classic() +
    labs(x = "ICC at household level",
       y = " GDP per capita (INR)",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(50000,100000,150000,200000,250000), limits=c(0, 300000)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/300000,expand=F)
stateobesehhfiggdp

###########figures for ICC on community level#####

#calculate p-values
summary(lm(GDP ~ ICC_diab_comm, data=stateprev.dat))
summary(lm(GDP ~ ICC_htn_comm, data=stateprev.dat))
summary(lm(GDP ~ ICC_csmoke_comm, data=stateprev.dat))
summary(lm(GDP ~ ICC_overweight_comm, data=stateprev.dat))
summary(lm(GDP ~ ICC_obese_comm, data=stateprev.dat))

statediabcommfiggdp <- stateprev.dat %>%
 
  ggplot(aes(y=GDP, x=ICC_diab_comm)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=GDP, x=ICC_diab_comm, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=GDP, x=ICC_diab_comm, label = state_lab)) +
 theme_classic() +
    labs(x = "ICC at community level",
       y = " GDP per capita (INR)",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(50000,100000,150000,200000,250000), limits=c(0, 300000)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/300000,expand=F)
statediabcommfiggdp

#######htn/comm

statehtncommfiggdp <- stateprev.dat %>%
 
  ggplot(aes(y=GDP, x=ICC_htn_comm)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=GDP, x=ICC_htn_comm, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=GDP, x=ICC_htn_comm, label = state_lab)) +
 theme_classic() +
    labs(x = "ICC at community level",
       y = " GDP per capita (INR)",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(50000,100000,150000,200000,250000), limits=c(0, 300000)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/300000,expand=F)
statehtncommfiggdp

###smoke/comm

statesmokecommfiggdp <- stateprev.dat %>%
 
  ggplot(aes(y=GDP, x=ICC_csmoke_comm)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=GDP, x=ICC_csmoke_comm, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=GDP, x=ICC_csmoke_comm, label = state_lab)) +
 theme_classic() +
   labs(x = "ICC at community level",
       y = " GDP per capita (INR)",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(50000,100000,150000,200000,250000), limits=c(0, 300000)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/300000,expand=F)
statesmokecommfiggdp

##overweight/comm

stateoverweightcommfiggdp <- stateprev.dat %>%
 
  ggplot(aes(y=GDP, x=ICC_overweight_comm)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=GDP, x=ICC_overweight_comm, color=zone, shape=zone), size=2.5) +
  geom_text_repel(aes(y=GDP, x=ICC_overweight_comm, label = state_lab)) +
  theme_classic() +
    labs(x = "ICC at community level",
       y = " GDP per capita (INR)",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(50000,100000,150000,200000,250000), limits=c(0, 300000)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/300000,expand=F)
stateoverweightcommfiggdp

##obese/comm

stateobesecommfiggdp <- stateprev.dat %>%
 
  ggplot(aes(y= GDP, x=ICC_obese_comm)) +
  geom_smooth(method='glm', se= FALSE, color="black") +
  geom_jitter(aes(y=GDP, x=ICC_obese_comm, color=zone, shape=zone), size=2.5) +
 geom_text_repel(aes(y=GDP, x=ICC_obese_comm, label = state_lab)) +
 theme_classic() +
  labs(x = "ICC at community level",
       y = " GDP per capita (INR)",
       fill="") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18, face="bold"),
        legend.text=element_text(size=18),
        legend.title = element_blank(),
        #legend.position="bottom",
        axis.title.x = element_text(margin = margin(t = 20)),
        axis.title.y = element_text(margin = margin(r = 20)),
        strip.text.x = element_text(size=20, face="bold"),
        strip.background = element_blank(),
        panel.spacing = unit(2, "lines")) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous(breaks = c(50000,100000,150000,200000,250000), limits=c(0, 300000)) +
  scale_x_continuous(breaks = c(0.1, 0.2, 0.3), limits=c(0, 0.4)) +
  coord_fixed(0.4/300000,expand=F)
stateobesecommfiggdp
```

```{r APPENDIX: non-fasted analysis}

#######CAUTION: now "india" will be the non-fasted version######

india<- filter(india.cvd, is.na(female)==F & is.na(bmi)==F & is.na(currsmoke)==F) #sex and BMI and csmoke
india <- 
  filter(india, is.na(htn_narrow_avg)==F&is.na(diab_narrow_nf)==F)

india$ex_htn_narrow_ind <-india$htn_narrow_avg
india$ex_diab_narrow_ind <- india$diab_narrow
india$csmoke <-india$currsmoke
india$wealth_quintile_rurb<-india$q_ai_nat_rurb
india$marital <- india$married
india$age_10yr <-india$age_grp

#################Fit the random effects model and calculate the ICC###############
#response must be numeric
#india$diab_narrow_nf<-as.numeric(levels(india$diab_narrow_nf))[india$diab_narrow_nf]
india$ex_htn_narrow_ind<-as.numeric(levels(india$ex_htn_narrow_ind))[india$ex_htn_narrow_ind]
india$csmoke<-as.numeric(levels(india$csmoke))[india$csmoke]
#india$overweight23<-as.numeric(levels(india$overweight23))[india$overweight23]
#india$obese27<-as.numeric(levels(india$obese27))[india$obese27]


#Table S4: Clustering of cardiovascular disease risk factors at the state, district, community, and household level assuming all AHS participants non-fasted

#state level################################################################
diab_state_nf <- lmer(diab_narrow_nf ~ 1 + (1|state), data=india)
calc.icc(diab_state_nf)

htn_state <- lmer(ex_htn_narrow_ind ~ 1 + (1|state), data=india)
calc.icc(htn_state) #has to be the same as with ex_diab_ind(same number of obs=1082100, so do not change)

csmoke_state <- lmer(csmoke ~ 1 + (1|state), data=india)
calc.icc(csmoke_state)#the same

overweight_state <- lmer(overweight23 ~ 1 + (1|state), data=india)
calc.icc(overweight_state)#the same

obese_state <- lmer(obese27 ~ 1 + (1|state), data=india)
calc.icc(obese_state)#the same


#calculcate confidence intervals via bootstrapping
#Calculate the bootstrap distribution
boot_diab_state_nf <- bootMer(diab_state_nf, calc.icc, nsim=500) #minimum simulations 500

#Draw from the bootstrap distribution the usual 95% upper and lower confidence limits
quantile(boot_diab_state_nf$t, c(0.025, 0.975))

#district level###############################################################
diab_dist_nf <- lmer(diab_narrow_nf ~ 1 + (1|state_dist), data=india)
calc.icc(diab_dist_nf)

#calculcate confidence intervals via bootstrapping
#Calculate the bootstrap distribution
memory.limit(size=50000)
boot_diab_dist_nf <- bootMer(diab_dist_nf, calc.icc, nsim=500) #minimum simulations 500


#Draw from the bootstrap distribution the usual 95% upper and lower confidence limits
quantile(boot_diab_dist_nf$t, c(0.025, 0.975))


##########################PSU level###############################################
diab_psu_nf <- lmer(diab_narrow_nf ~ 1 + (1|psuid), data=india)
calc.icc(diab_psu_nf)


#calculcate confidence intervals via bootstrapping
#Calculate the bootstrap distribution
boot_diab_psu_nf <- bootMer(diab_psu_nf, calc.icc, nsim=500)
#minimum simulations 500

#Draw from the bootstrap distribution the usual 95% upper and lower confidence limits
quantile(boot_diab_psu_nf$t, c(0.025, 0.975))


####household level###################################################################
diab_hh_nf <- lmer(diab_narrow_nf ~ 1 + (1|hhid), data=india)
calc.icc(diab_hh_nf)

#calculcate confidence intervals via bootstrapping
#Calculate the bootstrap distribution
boot_diab_hh_nf <- bootMer(diab_hh_nf, calc.icc, nsim=500) #minimum simulations 500

#Draw from the bootstrap distribution the usual 95% upper and lower confidence limits
quantile(boot_diab_hh_nf$t, c(0.025, 0.975))

##############Stratifying by state#############

india$ex_diab_narrow_ind_ahs_notfasted <- india$diab_narrow_nf
india$hh_id <-india$hhid
india$hh_id<-as.numeric(levels(india$hh_id))[india$hh_id]

#divide dataset into dataset by state
AndamanandNicobarIslands <- india %>%
  filter(state== "Andaman and Nicobar")
AndhraPradesh <- india %>%
  filter(state== "Andhra Pradesh")
ArunachalPradesh <- india %>%
  filter(state== "Arunachal Pradesh")
Assam <- india %>%
  filter(state== "Assam")
Bihar <- india %>%
  filter(state== "Bihar")
Chandigarh <- india %>%
  filter(state== "Chandigarh")
Chhattisgarh <- india %>%
  filter(state== "Chhattisgarh")
#DadraandNagarHaveli <- india %>%
  #filter(state== "Dadra and Nagar Haveli")
DamanandDiu <- india %>%
  filter(state== "Daman and Diu")
Delhi <- india %>%
  filter(state== "Delhi")
Goa <- india %>%
  filter(state== "Goa")
#Gujarat <- india %>%
  #filter(state== "Gujarat")
Haryana <- india %>%
  filter(state== "Haryana")
HimachalPradesh <- india %>%
  filter(state== "Himachal Pradesh")
#JammuandKashmir <- india %>%
  #filter(state== "Jammu and Kashmir")
Jharkhand <- india %>%
  filter(state== "Jharkhand")
Karnataka <- india %>%
  filter(state== "Karnataka")
Kerala <- india %>%
  filter(state== "Kerala")
#Lakshadweep <- india %>%
  #filter(state== "Lakshadweep")
MadhyaPradesh <- india %>%
  filter(state== "Madhya Pradesh")
Maharashtra<- india %>%
  filter(state== "Maharashtra")
Manipur <- india %>%
  filter(state== "Manipur")
Meghalaya <- india %>%
  filter(state== "Meghalaya")
Mizoram <- india %>%
  filter(state== "Mizoram")
Nagaland <- india %>%
  filter(state== "Nagaland")
Odisha <- india %>%
  filter(state== "Odisha")
Puducherry <- india %>%
  filter(state== "Puducherry")
Punjab <- india %>%
  filter(state== "Punjab")
Rajasthan <- india %>%
  filter(state== "Rajasthan")
Sikkim <- india %>%
  filter(state== "Sikkim")
TamilNadu <- india %>%
  filter(state== "Tamil Nadu")
Telangana <- india %>%
  filter(state== "Telangana")
Tripura <- india %>%
  filter(state== "Tripura")
UttarPradesh <- india %>%
  filter(state== "Uttar Pradesh")
Uttarakhand<- india %>%
  filter(state== "Uttarakhand")
WestBengal <- india %>%
  filter(state== "West Bengal")

#Table S5: Intracluster correlation coefficients at the household level by state assuming all AHS participants non-fasted	

#DIABETES 
diab_hh_nf_aani <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data= AndamanandNicobarIslands)
calc.icc(diab_hh_nf_aani)
diab_hh_nf_ap <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=AndhraPradesh)
calc.icc(diab_hh_nf_ap)
diab_hh_nf_arp <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=ArunachalPradesh)
calc.icc(diab_hh_nf_arp)
diab_hh_nf_as <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Assam)
calc.icc(diab_hh_nf_as)
diab_hh_nf_bi <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Bihar)
calc.icc(diab_hh_nf_bi)
diab_hh_nf_ch <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Chandigarh)
calc.icc(diab_hh_nf_ch)
diab_hh_nf_chh <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Chhattisgarh)
calc.icc(diab_hh_nf_chh)
#diab_hh_nf_danh <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=DadraandNagarHaveli)
#calc.icc(diab_hh_nf_danh)
diab_hh_nf_dad <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=DamanandDiu)
calc.icc(diab_hh_nf_dad)
diab_hh_nf_de <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Delhi)
calc.icc(diab_hh_nf_de)
diab_hh_nf_goa <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Goa)
calc.icc(diab_hh_nf_goa)
#diab_hh_nf_gu <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Gujarat)
#calc.icc(diab_hh_nf_gu)
diab_hh_nf_ha <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Haryana)
calc.icc(diab_hh_nf_ha)
diab_hh_nf_hp <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=HimachalPradesh)
calc.icc(diab_hh_nf_hp)
#diab_hh_nf_jk <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=JammuandKashmir)
#calc.icc(diab_hh_nf_jk)
diab_hh_nf_jh <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Jharkhand)
calc.icc(diab_hh_nf_jh)
diab_hh_nf_ka <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Karnataka)
calc.icc(diab_hh_nf_ka)
diab_hh_nf_ke <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Kerala)
calc.icc(diab_hh_nf_ke)
#diab_hh_nf_la <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Lakshadweep)
#calc.icc(diab_hh_nf_la)
diab_hh_nf_mp <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=MadhyaPradesh)
calc.icc(diab_hh_nf_mp)
diab_hh_nf_mah <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Maharashtra)
calc.icc(diab_hh_nf_mah)
diab_hh_nf_man <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Manipur)
calc.icc(diab_hh_nf_man)
diab_hh_nf_me <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Meghalaya)
calc.icc(diab_hh_nf_me)
diab_hh_nf_mi <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Mizoram)
calc.icc(diab_hh_nf_mi)
diab_hh_nf_na <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Nagaland)
calc.icc(diab_hh_nf_na)
diab_hh_nf_od <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Odisha)
calc.icc(diab_hh_nf_od)
diab_hh_nf_pud <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Puducherry)
calc.icc(diab_hh_nf_pud)
diab_hh_nf_pun <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Punjab)
calc.icc(diab_hh_nf_pun)
diab_hh_nf_ra <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Rajasthan)
calc.icc(diab_hh_nf_ra)
diab_hh_nf_si <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Sikkim)
calc.icc(diab_hh_nf_si)
diab_hh_nf_tn <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=TamilNadu)
calc.icc(diab_hh_nf_tn)
diab_hh_nf_te <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Telangana)
calc.icc(diab_hh_nf_te)
diab_hh_nf_tri <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Tripura)
calc.icc(diab_hh_nf_tri)
diab_hh_nf_up <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=UttarPradesh)
calc.icc(diab_hh_nf_up)
diab_hh_nf_ut <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=Uttarakhand)
calc.icc(diab_hh_nf_ut)
diab_hh_nf_wb <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|hh_id), data=WestBengal)
calc.icc(diab_hh_nf_wb)

############calculate bootstraps for state iccs########

boot_diab_hh_nf_aani <- bootMer(diab_hh_nf_aani, calc.icc, nsim=500)
boot_diab_hh_nf_ap <- bootMer(diab_hh_nf_ap, calc.icc, nsim=500)
boot_diab_hh_nf_arp <- bootMer(diab_hh_nf_arp, calc.icc, nsim=500)
boot_diab_hh_nf_as <- bootMer(diab_hh_nf_as, calc.icc, nsim=500)
boot_diab_hh_nf_bi <- bootMer(diab_hh_nf_bi, calc.icc, nsim=500)
boot_diab_hh_nf_ch <- bootMer(diab_hh_nf_ch, calc.icc, nsim=500)
boot_diab_hh_nf_chh <- bootMer(diab_hh_nf_chh, calc.icc, nsim=500)
#boot_diab_hh_nf_danh <- bootMer(diab_hh_nf_danh, calc.icc, nsim=500)
boot_diab_hh_nf_dad <- bootMer(diab_hh_nf_dad, calc.icc, nsim=500)
boot_diab_hh_nf_de <- bootMer(diab_hh_nf_de, calc.icc, nsim=500) 
boot_diab_hh_nf_goa <- bootMer(diab_hh_nf_goa, calc.icc, nsim=500)
#boot_diab_hh_nf_gu <- bootMer(diab_hh_nf_gu, calc.icc, nsim=500)
boot_diab_hh_nf_ha <- bootMer(diab_hh_nf_ha, calc.icc, nsim=500)
boot_diab_hh_nf_hp <- bootMer(diab_hh_nf_hp, calc.icc, nsim=500)
#boot_diab_hh_nf_jk <- bootMer(diab_hh_nf_jk, calc.icc, nsim=500)
boot_diab_hh_nf_jh <- bootMer(diab_hh_nf_jh, calc.icc, nsim=500)
boot_diab_hh_nf_ka <- bootMer(diab_hh_nf_ka, calc.icc, nsim=500)
boot_diab_hh_nf_ke <- bootMer(diab_hh_nf_ke, calc.icc, nsim=500)
#boot_diab_hh_nf_la <- bootMer(diab_hh_nf_la, calc.icc, nsim=500)
boot_diab_hh_nf_mp <- bootMer(diab_hh_nf_mp, calc.icc, nsim=500)
boot_diab_hh_nf_mah <- bootMer(diab_hh_nf_mah, calc.icc, nsim=500)
boot_diab_hh_nf_man <- bootMer(diab_hh_nf_man, calc.icc, nsim=500)
boot_diab_hh_nf_me <- bootMer(diab_hh_nf_me, calc.icc, nsim=500)
boot_diab_hh_nf_mi <- bootMer(diab_hh_nf_mi, calc.icc, nsim=500) 
boot_diab_hh_nf_na <- bootMer(diab_hh_nf_na, calc.icc, nsim=500)
boot_diab_hh_nf_od <- bootMer(diab_hh_nf_od, calc.icc, nsim=500)
boot_diab_hh_nf_pud <- bootMer(diab_hh_nf_pud, calc.icc, nsim=500)
boot_diab_hh_nf_pun <- bootMer(diab_hh_nf_pun, calc.icc, nsim=500)
boot_diab_hh_nf_ra <- bootMer(diab_hh_nf_ra, calc.icc, nsim=500)
boot_diab_hh_nf_si <- bootMer(diab_hh_nf_si, calc.icc, nsim=500)
boot_diab_hh_nf_tn <- bootMer(diab_hh_nf_tn, calc.icc, nsim=500)
boot_diab_hh_nf_te <- bootMer(diab_hh_nf_te, calc.icc, nsim=500)
boot_diab_hh_nf_tri <- bootMer(diab_hh_nf_tri, calc.icc, nsim=500)
boot_diab_hh_nf_up <- bootMer(diab_hh_nf_up, calc.icc, nsim=500)
boot_diab_hh_nf_ut <- bootMer(diab_hh_nf_ut, calc.icc, nsim=500)
boot_diab_hh_nf_wb <- bootMer(diab_hh_nf_wb, calc.icc, nsim=500)

############################################

quantile(boot_diab_hh_nf_aani$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_ap$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_arp$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_as$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_bi$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_ch$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_chh$t, c(0.025, 0.975)) 
#quantile(boot_diab_hh_nf_danh$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_dad$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_de$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_goa$t, c(0.025, 0.975)) 
#quantile(boot_diab_hh_nf_gu$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_ha$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_hp$t, c(0.025, 0.975)) 
#quantile(boot_diab_hh_nf_jk$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_jh$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_ka$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_ke$t, c(0.025, 0.975)) 
#quantile(boot_diab_hh_nf_la$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_mp$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_mah$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_man$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_me$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_mi$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_na$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_od$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_pud$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_pun$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_ra$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_si$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_tn$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_te$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_tri$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_up$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_ut$t, c(0.025, 0.975)) 
quantile(boot_diab_hh_nf_wb$t, c(0.025, 0.975))

##############community level#################################################

#Table S6: Intracluster correlation coefficients at the community level by state assuming all AHS participants non-fasted

diab_psu_nf_aani <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data= AndamanandNicobarIslands)
calc.icc(diab_psu_nf_aani)
diab_psu_nf_ap <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=AndhraPradesh)
calc.icc(diab_psu_nf_ap)
diab_psu_nf_arp <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=ArunachalPradesh)
calc.icc(diab_psu_nf_arp)
diab_psu_nf_as <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Assam)
calc.icc(diab_psu_nf_as)
diab_psu_nf_bi <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Bihar)
calc.icc(diab_psu_nf_bi)
diab_psu_nf_ch <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Chandigarh)
calc.icc(diab_psu_nf_ch)
diab_psu_nf_chh <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Chhattisgarh)
calc.icc(diab_psu_nf_chh)
#diab_psu_nf_danh <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=DadraandNagarHaveli)
#calc.icc(diab_psu_nf_danh)
diab_psu_nf_dad <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=DamanandDiu)
calc.icc(diab_psu_nf_dad)
diab_psu_nf_de <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Delhi)
calc.icc(diab_psu_nf_de)
diab_psu_nf_goa <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Goa)
calc.icc(diab_psu_nf_goa)
#diab_psu_nf_gu <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Gujarat)
#calc.icc(diab_psu_nf_gu)
diab_psu_nf_ha <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Haryana)
calc.icc(diab_psu_nf_ha)
diab_psu_nf_hp <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=HimachalPradesh)
calc.icc(diab_psu_nf_hp)
#diab_psu_nf_jk <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=JammuandKashmir)
#calc.icc(diab_psu_nf_jk)
diab_psu_nf_jh <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Jharkhand)
calc.icc(diab_psu_nf_jh)
diab_psu_nf_ka <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Karnataka)
calc.icc(diab_psu_nf_ka)
diab_psu_nf_ke <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Kerala)
calc.icc(diab_psu_nf_ke)
#diab_psu_nf_la <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Lakshadweep)
#calc.icc(diab_psu_nf_la)
diab_psu_nf_mp <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=MadhyaPradesh)
calc.icc(diab_psu_nf_mp)
diab_psu_nf_mah <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Maharashtra)
calc.icc(diab_psu_nf_mah)
diab_psu_nf_man <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Manipur)
calc.icc(diab_psu_nf_man)
diab_psu_nf_me <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Meghalaya)
calc.icc(diab_psu_nf_me)
diab_psu_nf_mi <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Mizoram)
calc.icc(diab_psu_nf_mi)
diab_psu_nf_na <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Nagaland)
calc.icc(diab_psu_nf_na)
diab_psu_nf_od <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Odisha)
calc.icc(diab_psu_nf_od)
diab_psu_nf_pud <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Puducherry)
calc.icc(diab_psu_nf_pud)
diab_psu_nf_pun <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Punjab)
calc.icc(diab_psu_nf_pun)
diab_psu_nf_ra <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Rajasthan)
calc.icc(diab_psu_nf_ra)
diab_psu_nf_si <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Sikkim)
calc.icc(diab_psu_nf_si)
diab_psu_nf_tn <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=TamilNadu)
calc.icc(diab_psu_nf_tn)
diab_psu_nf_te <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Telangana)
calc.icc(diab_psu_nf_te)
diab_psu_nf_tri <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Tripura)
calc.icc(diab_psu_nf_tri)
diab_psu_nf_up <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=UttarPradesh)
calc.icc(diab_psu_nf_up)
diab_psu_nf_ut <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=Uttarakhand)
calc.icc(diab_psu_nf_ut)
diab_psu_nf_wb <- lmer(ex_diab_narrow_ind_ahs_notfasted ~ 1 + (1|psuid), data=WestBengal)
calc.icc(diab_psu_nf_wb)




############calculate bootstraps for state iccs for community clustering########

boot_diab_psu_nf_aani <- bootMer(diab_psu_nf_aani, calc.icc, nsim=500)
boot_diab_psu_nf_ap <- bootMer(diab_psu_nf_ap, calc.icc, nsim=500)
boot_diab_psu_nf_arp <- bootMer(diab_psu_nf_arp, calc.icc, nsim=500)
boot_diab_psu_nf_as <- bootMer(diab_psu_nf_as, calc.icc, nsim=500)
boot_diab_psu_nf_bi <- bootMer(diab_psu_nf_bi, calc.icc, nsim=500)
boot_diab_psu_nf_ch <- bootMer(diab_psu_nf_ch, calc.icc, nsim=500)
boot_diab_psu_nf_chh <- bootMer(diab_psu_nf_chh, calc.icc, nsim=500)
#boot_diab_psu_nf_danh <- bootMer(diab_psu_nf_danh, calc.icc, nsim=500)
boot_diab_psu_nf_dad <- bootMer(diab_psu_nf_dad, calc.icc, nsim=500)
boot_diab_psu_nf_de <- bootMer(diab_psu_nf_de, calc.icc, nsim=500) 
boot_diab_psu_nf_goa <- bootMer(diab_psu_nf_goa, calc.icc, nsim=500)
#boot_diab_psu_nf_gu <- bootMer(diab_psu_nf_gu, calc.icc, nsim=500)
boot_diab_psu_nf_ha <- bootMer(diab_psu_nf_ha, calc.icc, nsim=500)
boot_diab_psu_nf_hp <- bootMer(diab_psu_nf_hp, calc.icc, nsim=500)
#boot_diab_psu_nf_jk <- bootMer(diab_psu_nf_jk, calc.icc, nsim=500)
boot_diab_psu_nf_jh <- bootMer(diab_psu_nf_jh, calc.icc, nsim=500)
boot_diab_psu_nf_ka <- bootMer(diab_psu_nf_ka, calc.icc, nsim=500)
boot_diab_psu_nf_ke <- bootMer(diab_psu_nf_ke, calc.icc, nsim=500)
#boot_diab_psu_nf_la <- bootMer(diab_psu_nf_la, calc.icc, nsim=500)
boot_diab_psu_nf_mp <- bootMer(diab_psu_nf_mp, calc.icc, nsim=500)
boot_diab_psu_nf_mah <- bootMer(diab_psu_nf_mah, calc.icc, nsim=500)
boot_diab_psu_nf_man <- bootMer(diab_psu_nf_man, calc.icc, nsim=500)
boot_diab_psu_nf_me <- bootMer(diab_psu_nf_me, calc.icc, nsim=500)
boot_diab_psu_nf_mi <- bootMer(diab_psu_nf_mi, calc.icc, nsim=500) 
boot_diab_psu_nf_na <- bootMer(diab_psu_nf_na, calc.icc, nsim=500)
boot_diab_psu_nf_od <- bootMer(diab_psu_nf_od, calc.icc, nsim=500)
boot_diab_psu_nf_pud <- bootMer(diab_psu_nf_pud, calc.icc, nsim=500)
boot_diab_psu_nf_pun <- bootMer(diab_psu_nf_pun, calc.icc, nsim=500)
boot_diab_psu_nf_ra <- bootMer(diab_psu_nf_ra, calc.icc, nsim=500)
boot_diab_psu_nf_si <- bootMer(diab_psu_nf_si, calc.icc, nsim=500)
boot_diab_psu_nf_tn <- bootMer(diab_psu_nf_tn, calc.icc, nsim=500)
boot_diab_psu_nf_te <- bootMer(diab_psu_nf_te, calc.icc, nsim=500)
boot_diab_psu_nf_tri <- bootMer(diab_psu_nf_tri, calc.icc, nsim=500)
boot_diab_psu_nf_up <- bootMer(diab_psu_nf_up, calc.icc, nsim=500)
boot_diab_psu_nf_ut <- bootMer(diab_psu_nf_ut, calc.icc, nsim=500)
boot_diab_psu_nf_wb <- bootMer(diab_psu_nf_wb, calc.icc, nsim=500)


############################################

quantile(boot_diab_psu_nf_aani$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_ap$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_arp$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_as$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_bi$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_ch$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_chh$t, c(0.025, 0.975)) 
#quantile(boot_diab_psu_nf_danh$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_dad$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_de$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_goa$t, c(0.025, 0.975)) 
#quantile(boot_diab_psu_nf_gu$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_ha$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_hp$t, c(0.025, 0.975)) 
#quantile(boot_diab_psu_nf_jk$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_jh$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_ka$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_ke$t, c(0.025, 0.975)) 
#quantile(boot_diab_psu_nf_la$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_mp$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_mah$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_man$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_me$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_mi$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_na$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_od$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_pud$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_pun$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_ra$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_si$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_tn$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_te$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_tri$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_up$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_ut$t, c(0.025, 0.975)) 
quantile(boot_diab_psu_nf_wb$t, c(0.025, 0.975))

	

```


